Alexandre Puyguiraud

Philippe Gouze

Marco Dentz
email: marco.dentz@csic.es

Temporal evolution of solute dispersion in three-dimensional porous rocks

Keywords:

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-like medium, 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. Solute dispersion is quantified by effective and ensemble dispersion coefficients. The former is a measure for the typical width of the plume, the latter for the deformation, that is, the spread of the mixing front. Both dispersion coefficients evolve from the molecular diffusion coefficients toward a common finite asymptotic value. Their evolution is governed by the interplay between diffusion, pore-scale velocity fluctuations and the medium structure, which determine the characteristic diffusion and advection time scales. Dispersion in the sand-like medium evolves on the transverse diffusion time across a characteristic streamtube diameter, which is the mechanism by which pore-scale flow variability is sampled by the solute. Dispersion in the Berea sandstone in contrast is governed by both the diffusion time across a typical streamtube, and the diffusion time along a pore conduit. These insights shed light on the evolution of mixing fronts in porous rocks, with implications for the understanding and modeling of transport phenomena of microbes and reactive solutes in porous media.

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 (Datta et al., 2013;Alim et al., 2017;[START_REF] Valocchi | Mixing-limited reactions in porous media[END_REF]Puyguiraud et al., 2021) while the later is ubiquitously controlled by the concentration gradients. The heterogeneity of the porous medium that triggers the flow speed distribution is therefore a primary parameter that controls dispersion from pre-asymptotic to Fickian regime (Dentz et al., 2004;Sherman et al., 2021). Transport in porous media is considered in many fields of academic and industrial applications from materials 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 (Saffman, 1959;[START_REF] Whitaker | Diffusion and dispersion in porous media[END_REF]Gelhar & Axness, 1983;Dagan, 1990;Dentz et al., 2023). Here we focus on hydrodynamic dispersion due 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 (Brenner, 1980;Ahmadi et al., 1998;Koch & Brady, 1985;Scheven, 2013;Bijeljic & Blunt, 2006;Le Borgne et al., 2011;Souzy et al., 2020;Puyguiraud et al., 2021). Similar challenges are encountered for reactive transport that is controlled by the time resolved distribution of the solutes and their mixing. If the reaction thermodynamics and kinetics are known, then the goal is to be able to model the local reaction rate from knowing dispersion properties (Battiato et al., 2009;Battiato & Tartakovsky, 2011). However, it is well known that the advection-dispersion equation parameterized by constant asymptotic dispersion coefficients are not suited to evaluate the effective reaction rates, because it assumes full mixing whereas incomplete mixing is the rule during the pre-asymptotic (non-Fickian) dispersion regimes (Rolle et al., 2009;Le Borgne et al., 2010;Dentz et al., 2011;Le Borgne et al., 2011;Puyguiraud et al., 2021). Nevertheless, diffusion and transverse mixing tend to homogenize concentration and full mixing can be expected in the asymptotic regime, as long as the characteristic length of heterogeneity is finite. Clearly, the convergence rate toward asymptotic dispersion and full mixing depend on the medium heterogeneity, but characterizing the relationship is still challenging and requires investigating both mixing and spreading mechanisms at all scales.

Solute 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 (Hulin & Plona, 1989;Khrapitchev & Callaghan, 2003;Bijeljic et al., 2004;Gouze et al., 2021;Puyguiraud et al., 2021;Gouze et al., 2023). These studies, accounting for the heterogeneity as a whole, show that the pore structure shapes the evolution of dispersion during the pre-asymptotic regime and then determine 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-asymtotic 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 or line source, is called ensemble dispersion in the following. As outlined above, at preasymptotic times ensemble dispersion measures advective solute spreading rather than mixing. In fact, it measures the center of mass fluctuations of the partial plume that evolves from the point injections that constitute the spatially extended initial distribution (Bouchaud & Georges, 1990). Several authors studied these dispersion concept in the context of mixing and dispersion in porous media on the continuum scale characterized by spatially variable hydraulic conductivity (Attinger et al., 1999;Dentz et al., 2000;Fiori, 2001;Fiori & Dagan, 2000;[START_REF] Vanderborght | Concentration variance and spatial covariance in secondorder stationary heterogeneous conductivity fields[END_REF]Dentz & de Barros, 2015;De Barros et al., 2015;de Barros & Dentz, 2016). Dentz et al. (2000) analyzed the time evolution of the effective and ensemble dispersion coefficients. They showed that the time resolved ensemble dispersion coefficient is usually larger than the effective dispersion until the effective dispersion growth rate increases due transverse local dispersion and diffusion and eventually converges with the ensemble dispersion coefficient. This increase of the effective dispersion value denotes the convergence of average local mixing toward macroscopic mixing that accounts for heterogeneity as a whole. Because it is a quantitative way to discriminate mixing from spreading, the notion of effective dispersion has been discussed and used by several authors for the modeling of experimental and numerical reactive transport data (Cirpka, 2002;Jose et al., 2004;Perez et al., 2019Perez et al., , 2020;;Puyguiraud et al., 2020).

As discussed above, most works that analyze effective and ensemble dispersion to quan-tify the impact of spatial heterogeneity on solute mixing and spreading consider continuum scale fluctuations of the hydraulic conductivity. To the best of our knowledge, the concept of effective dispersion has not been studied for transport in three-dimensional porous media despite its potential to explain the overestimation of pore-scale mixing and reaction by constant asymptotic hydrodynamic dispersion coefficients (Kapoor et al., 1998;Gramling et al., 2002;Perez et al., 2019).

In the present communication we investigate in detail the temporal evolution of mixing and spreading mechanisms occurring in porous media, in order to evaluate the different regimes in relation with the porous media structure. To this end, we perform threedimensional direct numerical simulations of pore-scale flow and solute transport in a sandpack medium and in a Berea sandstone of distinctly different heterogeneity levels, that can be measured, for instance, by the variance the logarithm of the flow velocity distribution. Solute dispersion is quantified by the temporal evolution of the effective and of the ensemble dispersion coefficients. This paper is organized as follows: the methodology used to calculate flow and transport and measure dispersion are presented in Section 2. In Section 3, we present the analyze of the dispersion behavior in the sand pack and Berea samples and discuss how these information can help us depicting the different dispersion stages in relation with the porous media structure. Section 4 presents the conclusions of the study.

Methodology

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

∇ 2 u(x) = - 1 µ ∇p(x), ∇ • u(x) = 0, ( 1
)
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

+ ∇ • [u(x) -D∇] c(x, t) = 0, (2)
where c(x, t) is the solute concentration at position x and time t, and D is the molecular diffusion coefficient. The advection-diffusion equation ( 2) is equivalent to the Langevin equation (Risken, 1996)

dx(t) dt = u[x(t)] + √ 2Dξ(t), (3)
where ξ(t) is a Gaussian white noise with mean ξ i = 0 and covariance ξ j (t)ξ k (t) = δ jk δ(tt ); δ jk is the Kronecker delta.

Mixing versus spreading

In this section, we discuss plume mixing versus spreading in terms of effective and ensemble dispersion coefficients. Then, we pose an approximate model based on the concept of effective dispersion to upscale pore-scale mixing to the continuum scale.

We analyze the mixing and dispersion of a solute by considering the concentration distribution c(x, t) for the normalized plane source

c(x, t = 0) = ρ(x) = φ -1 δ(x 1 ) I(x ∈ Ω f ) wh , (4)
where Ω f denotes the fluid domain and I(•) is the indicator function, which 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

Ω dxρ(x) = φ, (5)
where Ω denotes the bulk domain, that is, the union of fluid domain and solid domain.

The solute distribution can be decomposed into partial plumes g(x, t|x ) that satisfy Eq. ( 2)

for the initial conditions

g(x, t = 0|x ) = δ(x -x )I(x ∈ Ω f ). (6)
Then, we can write the concentration distribution c(x, t) as

c(x, t) = Ω dx ρ(x )g(x, t|x ). (7)
Note that g(x, t|y , z ) is the Green function of the transport problem. In the following, we define a surrogate model for the Green function using the concept of effective dispersion.

Effective and ensemble dispersion coefficients

In order to define effective and ensemble dispersion coefficients, we consider the moments of the Green function g(x, t|x ) and the concentration distribution c(x, t). The first and second moments of g(x, t|x ) are defined by

m i (t; x ) = dxx i g(x, t|x ), (8)
m ij (t; x ) = dxx i x j g(x, t|x ). ( 9
)
The first moments m i (t; x ) determine the center of mass position of g(x, t|x ). The second centered moments

κ ij (t; x ) = m (2) ij (t; x ) -m (1) i (t; x )m (1) j (t; x ) ( 10
)
are measures for the spatial extension of the Green function. The average of κ ij (t; x ) over all Green functions defines the effective second centered moment

κ eff ij (t) = dx ρ(x )κ ij (t; x ). ( 11
)
It is a measure for the average width of the Green function. The temporal rate of growth of κ eff ij (t) is given by the effective dispersion coefficients

D eff ij (t) = 1 2 d dt κ e ij (t), (12)
The effective dispersion coefficient measures the rate of growth of the spatial variance of a concentration distribution that evolves from a point-like initial condition.

In full analogy, we define the first and second moments of c(x, t) as

m i (t) = dxx i c(x, t) = dx ρ(x )m i (t; x ), (13)
m ij (t) = dxx i x j c(x, t) = dx ρ(x )m ij (t; x ). ( 14
)
As per the second equality signs, the moments are determined by taking ensemble averages over the moments of the set of Green functions and as such are named the ensemble moments in the following. The second centered ensemble moments are defined by

κ ens ij (t) = m ij (t) -m i (t)m j (t). ( 15
)
They are measures for the spatial extension of the concentration distribution, or equivalently for the ensemble of Green functions. The temporal rate of growth of the second centered ensemble moments is measured by the ensemble dispersion coefficients

D ens ij (t) = 1 2 d dt κ ens ij (t). ( 16
)
The difference between the ensemble and effective variances,

δκ m ij (t) = dx ρ(x ) m (1) i (t; x ) -m (1) i (t) m (1) j (t; x ) -m (1) j (t) , (17)
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 that constitute the solute plume

δD m ij (t) = 1 2 d dt δκ m ij (t). ( 18
)
In the following, we study the effective and ensemble dispersion coefficients as well as the center of mass fluctuations in streamwise direction, that is, for i = j = 1.

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 and fluid flow

We study two three-dimensional porous media of different complexity, (i) a Berea sandstone sample and (ii) a sand pack sample illustrated in Figure 1 The Berea sample displays a complex pore structure with a porosity of φ = 0.18, see also (Puyguiraud et al., 2021). This type of porous rock is considered to be a pertinent large-scale homogeneous proxy of high permeability sedimentary reservoirs (Churcher et al., 1991). The sand pack sample has a high porosity of φ = 0.37 with a more regular structure of the pore space. The sand-pack image (Sand Pack LV60C) was 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 difference between the two porous medium samples can be illustrated by the distribution of flow speeds (Alhashmi et al., 2016) shown in Figure 1. The flow heterogeneity is measured by the variance σ 2 f of the 10 -6

10 -4 10 -2 0 2 4 10 -7 10 -6 10 -5 10 -4 10 -3 10 -2 10 -1 10 0 10 1 pdf v/ v

Sand pack Berea

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

manuscript submitted to Water Resources Research natural logarithm f = ln v of the flow speed v. For the Berea sandstone sample, we obtain σ 2 f = 10, for the sand pack sample σ 2 f = 2, that is, the Berea sample is significantly more heterogeneous. The characteristic pore length scale is 0 = 1.5 × 10 -6 m both for the Berea and sand pack samples.

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. Firstly, it perfectly fits the voxels of the X-Ray tomography images, and secondly, 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.

Fluid flow in the pore space is solved numerically using the SIMPLE algorithm implemented in OpenFOAM [START_REF] Weller | A tensorial approach to computational continuum mechanics using object-oriented techniques[END_REF]. 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.

Once the solver has converged, the flow velocities are extracted at the centers of the interfaces of the mesh (that is, 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 P e = 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.

Random walk particle tracking

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

x(t + ∆t) = x(t) + u[x(t)]∆t + √ 2D∆tζ(t), (19)

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 P e = 2000. 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 . 18). 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

Center of mass

δD m 11 (t) = σ 2 0 t, (20)
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 11 (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 mechanisms 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 a characteristic flow rate, the decorrelation length becomes c = Q c / v . Thus, the time scale at which particles decorrelate is

τ c = 2 c D = Q c D 0 τ v . (21)
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 11 (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 11 (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.

The behavior in the second time regime is characterized by the transverse velocity sampling of particles that are initialized at moderate to high flow velocities on the one hand and the persistence of particles in low velocity conducts on the other hand, which gives rise to the observed sub-linear increase of δD m 11 (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 a low velocities along a conduct are cutoff 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

D ens 11 (t) = σ 2 0 t, (22)
where σ 2 0 is the initial velocity variance. It behaves in the same way as ∆D m 11 (t), see Eq. ( 20).

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 11 (t) first remain at the value of the molecular diffusion coefficient and then increase steeply due to shear dispersion.

Figures 4b and5b show that the increase of the effective dispersion coefficients occurs for high P e at earlier non-dimensional times than for low P e. This indicates that the shear rate does not scale linearly with u . In fact, a typical shear rate can be written as

γ = v γ , ( 23
)
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 . Thus, the characteristic shear time scale

τ γ = γ -1 ∝ τ v / v 1/2 .
This dependence explains the differences in the time behaviors of the effective dispersion coefficients for different P e.

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. 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 P e, 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 andd, they converge with D ens 11 (t).

Sand pack sample

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 11 (t)

and the convergence of D eff 11 (t) toward D ens 11 (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 11 (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 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 Eq. ( 21), which is proportional to τ v . This observation explains the temporal evolution of the ensemble and effective dispersion coefficients for t < τ v .

Berea sandstone sample

Conclusions

We investigate solute dispersion in three-dimensional porous rocks using detailed numerical simulations of pore-scale flow and transport. We consider a sand-like medium, We argue that it is first important to realize that solute dispersion evolves in time, and on time scales that are relevant for the understanding of transport phenomena of reactive solutes and microbes, for example. Second, it is important to realize that there is a conceptual and quantitative difference between solute spreading, as quantified by ensemble dispersion, and solute mixing, which is represented here by effective dispersion because it measures the typical rate of growth of the width of a partial plume that evolves from a point-like injection. The temporal evolution of effective dispersion from molecular diffusion to asymptotic hydrodynamic dispersion sheds light on the evolution of mixing fronts in porous media, and may explain phenomena of incomplete mixing observed for fast chemical reactions in porous media. 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 (Datta et al., 2013;Alim et al., 2017;[START_REF] Valocchi | Mixing-limited reactions in porous media[END_REF]Puyguiraud et al., 2021) while the later is ubiquitously controlled by the concentration gradients. The heterogeneity of the porous medium that triggers the flow speed distribution is therefore a primary parameter that controls dispersion from pre-asymptotic to Fickian regime (Dentz et al., 2004;Sherman et al., 2021). Transport in porous media is considered in many fields of academic and industrial applications from materials 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 (Saffman, 1959;[START_REF] Whitaker | Diffusion and dispersion in porous media[END_REF]Gelhar & Axness, 1983;Dagan, 1990;Dentz et al., 2023). Here we focus on hydrodynamic dispersion due 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 (Brenner, 1980;Ahmadi et al., 1998;Koch & Brady, 1985;Scheven, 2013;Bijeljic & Blunt, 2006;Le Borgne et al., 2011;Souzy et al., 2020;Puyguiraud et al., 2021). Similar challenges are encountered for reactive transport that is controlled by the time resolved distribution of the solutes and their mixing. If the reaction thermodynamics and kinetics are known, then the goal is to be able to model the local reaction rate from knowing dispersion properties (Battiato et al., 2009;Battiato & Tartakovsky, 2011). However, it is well known that the advection-dispersion equation parameterized by constant asymptotic dispersion coefficients are not suited to evaluate the effective reaction rates, because it assumes full mixing whereas incomplete mixing is the rule during the pre-asymptotic (non-Fickian) dispersion regimes (Rolle et al., 2009;Le Borgne et al., 2010;Dentz et al., 2011;Le Borgne et al., 2011;Puyguiraud et al., 2021). Nevertheless, diffusion and transverse mixing tend to homogenize concentration and full mixing can be expected in the asymptotic regime, as long as the characteristic length of heterogeneity is finite. Clearly, the convergence rate toward asymptotic dispersion and full mixing depend on the medium heterogeneity, but characterizing the relationship is still challenging and requires investigating both mixing and spreading mechanisms at all scales.

Solute 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 (Hulin & Plona, 1989;Khrapitchev & Callaghan, 2003;Bijeljic et al., 2004;Gouze et al., 2021;Puyguiraud et al., 2021;Gouze et al., 2023). These studies, accounting for the heterogeneity as a whole, show that the pore structure shapes the evolution of dispersion during the pre-asymptotic regime and then determine 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-asymtotic 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 or line source, is called ensemble dispersion in the following. As outlined above, at preasymptotic times ensemble dispersion measures advective solute spreading rather than mixing. In fact, it measures the center of mass fluctuations of the partial plume that evolves from the point injections that constitute the spatially extended initial distribution (Bouchaud & Georges, 1990). Several authors studied these dispersion concept in the context of mixing and dispersion in porous media on the continuum scale characterized by spatially variable hydraulic conductivity (Attinger et al., 1999;Dentz et al., 2000;Fiori, 2001;Fiori & Dagan, 2000;[START_REF] Vanderborght | Concentration variance and spatial covariance in secondorder stationary heterogeneous conductivity fields[END_REF]Dentz & de Barros, 2015;De Barros et al., 2015;de Barros & Dentz, 2016). Dentz et al. (2000) analyzed the time evolution of the effective and ensemble dispersion coefficients. They showed that the time resolved ensemble dispersion coefficient is usually larger than the effective dispersion until the effective dispersion growth rate increases due transverse local dispersion and diffusion and eventually converges with the ensemble dispersion coefficient. This increase of the effective dispersion value denotes the convergence of average local mixing toward macroscopic mixing that accounts for heterogeneity as a whole. Because it is a quantitative way to discriminate mixing from spreading, the notion of effective dispersion has been discussed and used by several authors for the modeling of experimental and numerical reactive transport data (Cirpka, 2002;Jose et al., 2004;Perez et al., 2019Perez et al., , 2020;;Puyguiraud et al., 2020).

As discussed above, most works that analyze effective and ensemble dispersion to quan-tify the impact of spatial heterogeneity on solute mixing and spreading consider continuum scale fluctuations of the hydraulic conductivity. To the best of our knowledge, the concept of effective dispersion has not been studied for transport in three-dimensional porous media despite its potential to explain the overestimation of pore-scale mixing and reaction by constant asymptotic hydrodynamic dispersion coefficients (Kapoor et al., 1998;Gramling et al., 2002;Perez et al., 2019).

In the present communication we investigate in detail the temporal evolution of mixing and spreading mechanisms occurring in porous media, in order to evaluate the different regimes in relation with the porous media structure. To this end, we perform three-

Methodology

2.1 Pore-scale flow and transport

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

∇ 2 u(x) = - 1 µ ∇p(x), ∇ • u(x) = 0, ( 1
)
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

+ ∇ • [u(x) -D∇] c(x, t) = 0, (2)
where c(x, t) is the solute concentration at position x and time t, and D is the molecular diffusion coefficient. The advection-diffusion equation ( 2) is equivalent to the Langevin equation (Risken, 1996)

dx(t) dt = u[x(t)] + √ 2Dξ(t), (3)
where ξ(t) is a Gaussian white noise with mean ξ i = 0 and covariance ξ j (t)ξ k (t) = δ jk δ(tt ); δ jk is the Kronecker delta.

Mixing versus spreading

In this section, we discuss plume mixing versus spreading in terms of effective and ensemble dispersion coefficients. Then, we pose an approximate model based on the concept of effective dispersion to upscale pore-scale mixing to the continuum scale.

We analyze the mixing and dispersion of a solute by considering the concentration distribution c(x, t) for the normalized plane source

c(x, t = 0) = ρ(x) = φ -1 δ(x 1 ) I(x ∈ Ω f ) wh , (4)
where Ω f denotes the fluid domain and I(•) is the indicator function, which 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

Ω dxρ(x) = φ, (5)
where Ω denotes the bulk domain, that is, the union of fluid domain and solid domain.

The solute distribution can be decomposed into partial plumes g(x, t|x ) that satisfy Eq. ( 2)

for the initial conditions

g(x, t = 0|x ) = δ(x -x )I(x ∈ Ω f ). ( 6
)
Then, we can write the concentration distribution c(x, t) as

c(x, t) = Ω dx ρ(x )g(x, t|x ). ( 7
)
Note that g(x, t|y , z ) is the Green function of the transport problem. In the following, we define a surrogate model for the Green function using the concept of effective dispersion.

Effective and ensemble dispersion coefficients

In order to define effective and ensemble dispersion coefficients, we consider the moments of the Green function g(x, t|x ) and the concentration distribution c(x, t). The first and second moments of g(x, t|x ) are defined by

m i (t; x ) = dxx i g(x, t|x ), (8)
m ij (t; x ) = dxx i x j g(x, t|x ). ( 9
)
The first moments m i (t; x ) determine the center of mass position of g(x, t|x ). The second centered moments

κ ij (t; x ) = m (2) ij (t; x ) -m (1) i (t; x )m (1) j (t; x ) ( 10
)
are measures for the spatial extension of the Green function. The average of κ ij (t; x ) over all Green functions defines the effective second centered moment

κ eff ij (t) = dx ρ(x )κ ij (t; x ). ( 11
)
It is a measure for the average width of the Green function. The temporal rate of growth of κ eff ij (t) is given by the effective dispersion coefficients

D eff ij (t) = 1 2 d dt κ e ij (t), (12)
The effective dispersion coefficient measures the rate of growth of the spatial variance of a concentration distribution that evolves from a point-like initial condition.

In full analogy, we define the first and second moments of c(x, t) as

m i (t) = dxx i c(x, t) = dx ρ(x )m i (t; x ), (13)
m ij (t) = dxx i x j c(x, t) = dx ρ(x )m ij (t; x ). ( 14
)
As per the second equality signs, the moments are determined by taking ensemble averages over the moments of the set of Green functions and as such are named the ensemble moments in the following. The second centered ensemble moments are defined by

κ ens ij (t) = m ij (t) -m i (t)m j (t). ( 15
)
They are measures for the spatial extension of the concentration distribution, or equivalently for the ensemble of Green functions. The temporal rate of growth of the second centered ensemble moments is measured by the ensemble dispersion coefficients

D ens ij (t) = 1 2 d dt κ ens ij (t). ( 16
)
The difference between the ensemble and effective variances,

δκ m ij (t) = dx ρ(x ) m (1) i (t; x ) -m (1) i (t) m
(1)

j (t; x ) -m (1) j (t) , (17)
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 that constitute the solute plume

δD m ij (t) = 1 2 d dt δκ m ij (t). (18)
In the following, we study the effective and ensemble dispersion coefficients as well as the center of mass fluctuations in streamwise direction, that is, for i = j = 1.

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 and fluid flow

We study two three-dimensional porous media of different complexity, (i) a Berea sandstone sample and (ii) a sand pack sample illustrated in Figure 1 The Berea sample displays a complex pore structure with a porosity of φ = 0.18, see also (Puyguiraud et al., 2021). This type of porous rock is considered to be a pertinent large-scale homogeneous proxy of high permeability sedimentary reservoirs (Churcher et al., 1991). The sand pack sample has a high porosity of φ = 0.37 with a more regular structure of the pore space. The sand-pack image (Sand Pack LV60C) was 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 difference between the two porous medium samples can be illustrated by the distribution of flow speeds (Alhashmi et al., 2016) shown in Figure 1. The flow heterogeneity is measured by the variance σ 2 f of the natural logarithm f = ln v of the flow speed v. For the Berea sandstone sample, we obtain σ 2 f = 10, for the sand pack sample σ 2 f = 2, that is, the Berea sample is significantly more heterogeneous. The characteristic pore length scale is 0 = 1.5 × 10 -6 m both for the Berea and sand pack samples.

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. Firstly, it perfectly fits the voxels of the X-Ray tomography images, and secondly, 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.

Fluid flow in the pore space is solved numerically using the SIMPLE algorithm implemented in OpenFOAM [START_REF] Weller | A tensorial approach to computational continuum mechanics using object-oriented techniques[END_REF]. 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.

Once the solver has converged, the flow velocities are extracted at the centers of the interfaces of the mesh (that is, 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 P e = 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.

Random walk particle tracking

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

x(t + ∆t) = x(t) + u[x(t)]∆t + √ 2D∆tζ(t), (19)

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

(t) = σ 2 0 t, (20)
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 11 (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 mechanisms 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 a characteristic flow rate, the decorrelation length becomes c = Q c / v . Thus, the time scale at which particles decorrelate is

τ c = 2 c D = Q c D 0 τ v . (21)
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 11 (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 11 (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.

The behavior in the second time regime is characterized by the transverse velocity sampling of particles that are initialized at moderate to high flow velocities on the one hand and the persistence of particles in low velocity conducts on the other hand, which gives rise to the observed sub-linear increase of δD m 11 (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 a low velocities along a conduct are cutoff 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

D ens 11 (t) = σ 2 0 t, (22)
where σ 2 0 is the initial velocity variance. It behaves in the same way as ∆D m 11 (t), see Eq. ( 20).

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 11 (t) first remain at the value of the molecular diffusion coefficient and then increase steeply due to shear dispersion.

Figures 4b and5b show that the increase of the effective dispersion coefficients occurs for high P e at earlier non-dimensional times than for low P e. This indicates that the shear rate does not scale linearly with u . In fact, a typical shear rate can be written as

γ = v γ , ( 23
)
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 . Thus, the characteristic shear time scale

τ γ = γ -1 ∝ τ v / v 1/2 .
This dependence explains the differences in the time behaviors of the effective dispersion coefficients for different P e.

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. 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 P e, 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 andd, they converge with D ens 11 (t).

Sand pack sample

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 11 (t)

and the convergence of D eff 11 (t) toward D ens 11 (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 11 (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 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 Eq. ( 21), which is proportional to τ v . This observation explains the temporal evolution of the ensemble and effective dispersion coefficients for t < τ v .

Berea sandstone sample

Conclusions

We investigate solute dispersion in three-dimensional porous rocks using detailed numerical simulations of pore-scale flow and transport. We consider a sand-like medium, ple than for the sand pack sample. Solute dispersion is quantified by effective and ensemble dispersion coefficients. The former is defined in terms of the spatial average of the second-centered moments of the partial plumes (Green functions) that constitute the global solute distribution. Ensemble dispersion coefficients are defined in terms of the second centered moments of the global solute plume. Thus, the effective dispersion coefficients can be seen as a measure for the typical width of a mixing front, while the ensemble dispersion coefficients are a measure for its deformation due to the flow variability within the initial plume. The mechanisms that cause hydrodynamic dispersion are pore-scale flow variability and molecular diffusion, and govern the evolution of both the effective and ensemble dispersion coefficients. They eventually converge toward the same asymptotic value, which quantifies the impact of spatial heterogeneity on large-scale mixing.

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 behave 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 assume their asymptotic long time value on the time scale for diffusion over a typical pore length. These behaviors can be traced back to the network-like medium structure in case of the Berea sample, and the strong connectivity of pores in the sand pack. Thus, the evolution of solute dispersion reflects the medium structure, which determines the microscopic mass transfer mechanisms. 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.

We argue that it is first important to realize that solute dispersion evolves in time, and on time scales that are relevant for the understanding of transport phenomena of reactive solutes and microbes, for example. Second, it is important to realize that there is a conceptual and quantitative difference between solute spreading, as quantified by ensemble dispersion, and solute mixing, which is represented here by effective dispersion because it measures the typical rate of growth of the width of a partial plume that evolves from a point-like injection. The temporal evolution of effective dispersion from molecular diffusion to asymptotic hydrodynamic dispersion sheds light on the evolution of mixing fronts in porous media, and may explain phenomena of incomplete mixing observed for fast chemical reactions in porous media.

  The average pore length 0 , the mean streamwise flow velocity 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 P e = τ D /τ v = v 0 /D.

  where ζ = (ζ 1 , ζ 2 , ζ 3 ). The ζ i are independent random variables that are uniformly distributed in [-√ 3, √ 3]. The central limit theorem ensures that the sum of these uniform random variables is Gaussian distributed with zero mean and unit variance. The particle velocities u[x(t)] are interpolated from the velocities at the voxel faces using the algorithm of Mostaghimi et al. (2012), which implements a quadratic interpolation in the void voxels that are in contact with the solid and thus guarantees an accurate representation of the flow field in the vicinity of the solid-void interface. The time step is variable and chosen such that the particle displacement at a given step is shorter than or equal to the side length of a voxel. The time step varies from ∆t = 10 -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, 1.5 × 10 7 particles are uniformly placed at a plane perpendicular to the mean flow direction, see Figure 2 for the Berea sandstone. A similar setup is used for the sand-pack. We consider this scenario for P e = 200 and P e = 2000.

Figure 2 .

 2 Figure 2. Snapshots of the conservative simulation for the Berea sandstone for P e = 2000 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

 3 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 11 (t) of the variance of the center of mass positions m 1 (t|x ) of partial plumes g(x, t|x ) defined by (18). The center of mass of the global plume moves with the mean

Figure 3 .

 3 Figure 3. Temporal evolution of the center of mass position of the (black solid line) global plume, and (orange dashed lines) selected partial plumes for the sand-pack with (top left) P e = 200 and (top right) P e = 2000, and the Berea sample with (bottom left) P e = 200 and (bottom right) P e = 2000. The dashed vertical lines denote (black) the advection time scale τv, (yellow and orange) the respective diffusion time scales τD.

Figures 4

 4 Figures4 and 5show 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 = P e -1 τ v , diffusion dominates over advection, and both the ensemble and effective dispersion coefficients are equal to the molecular diffusion coefficient D. 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

FiguresFigure 4 .

 4 Figures 4a-dshow 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 Figure4a.

Figures

  Figures 5a-d show the evolution of the ensemble and effective dispersion coefficientsfor the Berea sandstone sample. As seen in Figure5a, 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 cross-over toward the constant asymptotic long time values occurs on the diffusion time scale τ D . As discussed in Section 3.1.2, here the temporal decorrelation of low velocities is due to diffusion along pore channels with the characteristic time scale τ D(Puyguiraud et al., 2021). Similary, the convergence of the effective dispersion coefficient shown in Figure 5b occurs on the time scale τ D . The cross-over of the effective to the ensemble dispersion coefficients shown in Figures 5c and d occurs on the decorrelation time scale τ c , see Eq. (21). This time scale is set by transverse diffusion across streamtubes, which is the mechanisms by which particles that originate at the same initial position start decorrelating and sampling different flow velocities. The independent sampling of flow velocities along trajectories between different particles is the ensemble mechanism of dispersion as measured by the ensemble dispersion coefficients, and therefore effective and ensemble dispersion converge on the scale τ c .

Figure 5 .

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

S0169772220303235

  doi: https://doi.org/10.1016/j.jconhyd.2020.103734 Souzy, M., Lhuissier, H., Méheust, Y., Borgne, T. L., & Metzger, B. (2020, March). Velocity distributions, dispersion and stretching in three-dimensional porous media. Journal of Fluid Mechanics, 891 . Retrieved from https://doi.org/ 1 Introduction

  dimensional direct numerical simulations of pore-scale flow and solute transport in a sandpack medium and in a Berea sandstone of distinctly different heterogeneity levels, that can be measured, for instance, by the variance the logarithm of the flow velocity distribution. Solute dispersion is quantified by the temporal evolution of the effective and of the ensemble dispersion coefficients. This paper is organized as follows: the methodology used to calculate flow and transport and measure dispersion are presented in Section 2. In Section 3, we present the analyze of the dispersion behavior in the sand pack and Berea samples and discuss how these information can help us depicting the different dispersion stages in relation with the porous media structure. Section 4 presents the conclusions of the study.

  The average pore length 0 , the mean streamwise flow velocity 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 P e = τ D /τ v = v 0 /D.

Figure 1 .

 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 (5mm 3 ) and of (right) the Berea sandstone (1mm 3 ). The grey and blue colors represent the pore space and the solid phase, respectively.

  where ζ = (ζ 1 , ζ 2 , ζ 3 ). The ζ i are independent random variables that are uniformly distributed in [-√ 3, √ 3]. The central limit theorem ensures that the sum of these uniform random variables is Gaussian distributed with zero mean and unit variance. The particle velocities u[x(t)] are interpolated from the velocities at the voxel faces using the algorithm of Mostaghimi et al. (2012), which implements a quadratic interpolation in the void voxels that are in contact with the solid and thus guarantees an accurate representation of the flow field in the vicinity of the solid-void interface. The time step is variable and chosen such that the particle displacement at a given step is shorter than or equal to the side length of a voxel. The time step varies from ∆t = 10 -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, 1.5 × 10 7 particles are uniformly placed at a plane perpendicular to the mean flow direction, see Figure 2 for the Berea sandstone. A similar setup is used for the sand-pack. We consider this scenario for P e = 200 and P e = 2000.

Figure 2 .

 2 Figure 2. Snapshots of the conservative simulation for the Berea sandstone for P e = 2000 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 .

 3 Figure 3. Temporal evolution of the center of mass position of the (black solid line) global plume, and (orange dashed lines) selected partial plumes for the sand-pack with (top left) P e = 200 and (top right) P e = 2000, and the Berea sample with (bottom left) P e = 200 and (bottom right) P e = 2000. The dashed vertical lines denote (black) the advection time scale τv, (yellow and orange) the respective diffusion time scales τD.

Figures 4

 4 Figures4 and 5show 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 = P e -1 τ v , diffusion dominates over advection, and both the ensemble and effective dispersion coefficients are equal to the molecular diffusion coefficient D. 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

FiguresFigure 4 .

 4 Figures 4a-dshow 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 Figure4a.

Figures

  Figures 5a-d show the evolution of the ensemble and effective dispersion coefficientsfor the Berea sandstone sample. As seen in Figure5a, 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 cross-over toward the constant asymptotic long time values occurs on the diffusion time scale τ D . As discussed in Section 3.1.2, here the temporal decorrelation of low velocities is due to diffusion along pore channels with the characteristic time scale τ D(Puyguiraud et al., 2021). Similary, the convergence of the effective dispersion coefficient shown in Figure 5b occurs on the time scale τ D . The cross-over of the effective to the ensemble dispersion coefficients shown in Figures 5c and d occurs on the decorrelation time scale τ c , see Eq. (21). This time scale is set by transverse diffusion across streamtubes, which is the mechanisms by which particles that originate at the same initial position start decorrelating and sampling different flow velocities. The independent sampling of flow velocities along trajectories between different particles is the ensemble mechanism of dispersion as measured by the ensemble dispersion coefficients, and therefore effective and ensemble dispersion converge on the scale τ c .

Figure 5 .

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

S0169772220303235

  doi: https://doi.org/10.1016/j.jconhyd.2020.103734 Souzy, M., Lhuissier, H., Méheust, Y., Borgne, T. L., & Metzger, B. (2020, March). Velocity distributions, dispersion and stretching in three-dimensional porous media. Journal of Fluid Mechanics, 891 . Retrieved from https://doi.org/

Acknowledgments

This project has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Sk lodowska-Curie grant agreement No 899546. MD gratefully acknowledge the support of the Spanish Minstry of Science and Innovation through the project HydroPore (PID2019-106887GB-C31). The simulation data displayed in the figures can be downloaded at http://hdl.handle.net/10261/331188. Acknowledgments This project has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Sk lodowska-Curie grant agreement No 899546. MD gratefully acknowledge the support of the Spanish Minstry of Science and Innovation through the project HydroPore (PID2019-106887GB-C31). The simulation data displayed in the figures can be downloaded at http://hdl.handle.net/10261/331188.

Abstract

We study the temporal evolution of solute dispersion in three-dimensional porous rocks of different heterogeneity and pore structure. To this end, we perform direct numerical simulations of pore-scale flow and transport in a sand-like medium, 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. Solute dispersion is quantified by effective and ensemble dispersion coefficients. The former is a measure for the typical width of the plume, the latter for the deformation, that is, the spread of the mixing front. Both dispersion coefficients evolve from the molecular diffusion coefficients toward a common finite asymptotic value. Their evolution is governed by the interplay between diffusion, pore-scale velocity fluctuations and the medium structure, which determine the characteristic diffusion and advection time scales. Dispersion in the sand-like medium evolves on the transverse diffusion time across a characteristic streamtube diameter, which is the mechanism by which pore-scale flow variability is sampled by the solute. Dispersion in the Berea sandstone in contrast is governed by both the diffusion time across a typical streamtube, and the diffusion time along a pore conduit. These insights shed light on the evolution of mixing fronts in porous rocks, with implications for the understanding and modeling of transport phenomena of microbes and reactive solutes in porous media.