Impact of Synthetic Porous Medium Geometric Properties on Solute Transport Using Direct 3 D Pore-Scale Simulations

National Research Council of Italy, Water Research Institute, Area della Ricerca di Roma 1-Montelibretti, Strada Provinciale 35d, km 0.7 Montelibretti Roma, Italy Institute of Geochemistry and Petrology, ETH Zurich, Clausiusstrasse 25, CH-8092 Zurich, Switzerland Department of Earth, Environmental and Planetary Sciences, Brown University, Providence, 02912 RI, USA Department of Computational Hydrosystems Helmholtz, Centre for Environmental Research, UFZ Permoserstr. 15 04318 Leipzig, Germany


Introduction
A quantitative relation between the structure of porous media and its effect on solute transport is fundamental to our understanding of groundwater contamination dynamics and the development of pollutant remediation strategies [1].Transport processes in porous media have been traditionally studied through the parametrization of macroscale properties (e.g., hydraulic conductivity or hydrodynamic dispersion coefficients) by means of volume-averaging or upscaling methods over a representative elementary volume.
However, these upscaling approaches are not valid when considering length scales on the order of a single pore or the complexity of processes (i.e., physical, chemical, and biological) that take place at the pore scale [2].Developing averaging approaches from the pore-scale approach can provide further information on the macroscopic medium properties and then enhance the understanding of transport effects which exhibit at a greater scale (e.g., field or laboratory scale) [3][4][5][6].
The development of digital image processing led to a large expansion of the quantitative analyses of rock and soil structure.Various tomographic techniques are nowadays available for a three-dimensional non-invasive visualization of structural properties with a spatial resolution down to microns [7][8][9][10][11].Even with the availability of these techniques, the cost of imaging limits the number of samples that can be analyzed [12].As an alternative, a statistical approach to generate media with given imposed structural features can be a cost-effective means for producing large numbers of realizations that represent those features of natural samples deemed relevant.The two approaches that are commonly used to generate synthetic porous media can be divided into statistical models (e.g., multipoint statistics and autocorrelation functions) and object-based methods, where geometrical objects are randomly placed in a domain to simulate diagenesis processes [8].
The present research is aimed at analyzing the relationships between geometric features of different porous media and transport processes through direct simulations of fluid flow and advection-diffusion of non-reactive solute using the lattice Boltzmann method (LBM).The LBM is selected for its ability to easily model different phenomena and its suitability to be parallelized and to exploit the computing power of HPC systems [37].The question we address here is which geometric properties of the porous media control fluid flow and solute transport?In order to perform a sensitivity analysis on a wide class of porous microstructures, multiple synthetic porous media are generated using a geostatistical approach.We also consider two available sandpack samples (https://www.imperial.ac.uk/earth-science/research/research-groups/perm/research/ pore-scale-modelling/micro-ct-images-and-networks/, last access 16/01/2019) to compare results for synthetic and natural porous media.The choice of sandpack images is motivated because alluvial aquifers are the most commonly contaminated systems [1].Each porous medium is described by the Minkowski functionals: porosity, specific surface area, mean curvature, and Euler characteristic [38].

Materials and Methods
To identify the dominant pore-scale parameters that control the transport of non-reactive solutes, LBM simulations (Section 2.2) are run over synthetic porous media and micro-CT images of sandpacks (Section 2.1).A specific simulation setup is designed to first consider separately the impact of advective and diffusive transport (Section 2.3).Our objective is to reduce the number of variables to consider for the advection-diffusion simulations, which are subsequently carried out on a selection of porous media.In fact, a full sensitivity analysis between the correlation lengths (ξ) set for the generation model, the Minkowski functionals for porous media characterizations (porosity (φ), specific surface area (SSA), mean curvature (MC), and Euler number (χ V )), and tortuosity (τ) over different transport regimes would require a computational effort that goes beyond our current capacity.Finally, the metrics used to quantify the simulation results over the selected parameters are described in Section 2.4.
Acronyms and abbreviations adopted throughout the paper are summarized in Supplementary Material S3.
2.1.3D Porous Medium.The procedure to generate the porous media is described in details in Di Palma et al. [23].
Here, we only summarize the main steps as follows: (i) the 3D domain and resolution of the porous medium are set; (ii) the geostatistical model is assigned, setting the type of semivariogram (exponential and Gaussian) and the correlation length ξ.Pore-scale variability (i.e., for multiple length scale) is neglected; (iii) a 3D Gaussian random field (zero mean and variance of 1) following the assigned spatial autocorrelation function (semivariogram) is generated with the TFracGen code [39]; (iv) a threshold is imposed to the normal cumulated distribution function to obtain an indicator field of pores and grains with an assigned target porosity φ; and (v) the spatial autocorrelation of the generated indicator field is estimated by fitting the Matérn covariance model [40,41] as a function of the smoothness parameter ν m and the associated practical range.This step is necessary as the transformation of a Gaussian field into an indicator field does not preserve the characteristics of the spatial autocorrelation [12,42].The resulting smoothness parameters of the Matèrn model for the porous media generated by exponential and Gaussian variograms are ν m = 0 35 and ν m = 1 2, respectively.These values are consistent with the findings of Di Palma et al. [23].However, for the sake of simplicity, we refer for the rest of the paper to EXX and GXX to porous media generated by exponential and Gaussian functions, respectively, where XX refers to the correlation length of each medium.
In this work, each porous medium is generated at a spatial resolution of 10 μm for a domain size of 3 mm 3 , resulting in a number of voxels equal to 300 3 .A greater calculation domain would have involved a huge computationally cost, especially for advection-diffusion simulations.Porous media are generated in a range between 3% and 11% of the ratio correlation length over domain size to guarantee a suitable balance between the spatial resolution impact on short correlation lengths and the representativeness of porous medium generation for great correlation lengths compared to domain size [23].It is worth stressing that the synthetic porous media are isotropic (i.e., show only one independent length scale) and statistically homogeneous.We are aware that this modelling assumption constitutes a limitation of our approach, e.g., for sedimentary porous media that often show pronounced anisotropy between the horizontal and vertical directions.However, we think that this is justified by the reduction of 2 Geofluids the numbers of descriptors required to describe each medium.Such limitation does not significantly affect the interpretability of our results.In order to compare the analysis of the synthetic porous media with real samples, we performed additionally flow and solute transport simulations on two different sandpacks (i.e., F42C and LV60A) with the same spatial resolution of 10 μm and 300 3 voxels, whose micro-CT images are available from the pore-scale modelling Consortium of Imperial College (https://figshare.com/articles/F42C_sandstone/1189261,https://figshare.com/articles/LV60A_sandpack/1153795,last access 15/01/2019).Both the X-ray scanned sandpack media have been verified to be isotropic and homogenous [17].
To describe the spatial characteristics of the generated porous media, we adopted the four Minkowski functionals and the tortuosity.The Minkowski functionals provide a detailed description of the geometric characteristics of a porous medium and are widely used in mathematical morphology [12,38].We briefly summarize the general definitions for each of them but refer interested readers to Vogel et al. [38] for a complete treatment of the Minkowski functional.
The first functional is the porosity: defined as the volume of the voids V voids over the total volume of the medium (V).
The second functional is the specific surface area, as the interstitial surface area of the voids and solids: where A s is the surface area of the grains.The third functional is the mean curvature (MC) representing a measure of the average pore size of the sample over the total volume: Here, W 1 is related to the mean width and is defined as the mean value of the distance between a pair of parallel support planes [12,43].
Finally, the fourth functional is the total curvature, also known as the Euler number that measures the connectivity of the pore space.This is more easily defined by Vogel [38,44,45] as a topological parameter: where N are the isolated objects, C are interconnected regions or loops, and H are the completely enclosed cavities.Broadly speaking, a negative χ V corresponds to high pore space connectivity [45].We calculate the tortuosity following the definition of hydraulic tortuosity in Clennell [46] as τ = L e /L , L e being the average length of the fluid paths and L is the straightline length in the direction of flow [46].L e is obtained from steady-state flow simulations by counting the number of streamlines passing through the porous medium.Streamlines are calculated using MATLAB® function "stream3," which is based on an explicit Euler integration of the ordinary differential equation x t , y t , z t = u x t , y t , z t , v x t , y t , z t , w x t , y t , z t for a given set of starting points x t = 0 , y t = 0 , z t = 0 , , where u, v, and w represent the components of the velocity.
In our 3D media, the porosity φ and correlation length ξ are imposed in the generation process, whereas the Minkowski functionals and tortuosity are calculated after the porous medium realization and fluid flow simulations, respectively.
The imposed geometric properties of the generated porous media and associated nomenclature (ID) are reported in Table 1, together with the corresponding estimated geometric properties for the scanned samples.
Figure 1 shows the geometric descriptors of the generated and scanned porous media (both imposed or calculated using the Minkowski functionals).The dispersion of imposed properties is by construction regularly distributed over the imposed ranges of correlation length and porosity, including both micro-CT sandpacks (Figure 1(a)).The SSA values of the EXX are approximately two or three times greater than those of the GXX and almost independent of the porosity.Moreover, the SSA values of the micro-CT are closer to the GXX than to the EXX (Figure 1(b)).On the other hand, the SSA variability introduced by the investigated range of generated pore size is similar to that of EXX and GXX (Figure 1(c)).The dependence of the mean curvature (MC) on correlation length (Figure 1(d)) is similar to the dependence of the specific surface area (SSA) on the porosity (Figure 1(e)) while these two functionals appear correlated (Figure 1(f)).This implies that the two semivariogram models (Gaussian and exponential) adopted for synthetic generation lead to different average pore sizes for the same imposed correlation length.
The Euler number χ V values (Figures 1(g)-1(j)) are consistent with literature data for similar porous media [38], showing overall negative values, except for one sample.A good correlation is found between χ V and SSA (Figure 1(i)), similar to the second versus the third Minkowski functional (Figure 1(f)).Indeed, the Euler number and the mean curvature seem to be correlated (Figure 1(j)).The highest values of the χ V for the media generated using the exponential semivariogram model and the lowest correlation lengths (e.g., upper left black triangles in Figure 1(g)) are due to the presence of many isolated small structures, which result in a massive inflation of the connectivity.The issue is an artifact of the binarization procedure, which can be avoided by using a finer spatial resolution for the exponential semivariogram model or a dedicated postprocessing routine.
A principal component analysis has been carried out on the set of descriptors presented in Figure 1 (not shown here).The loads of the first PCA component confirm the strong relationship among SSA, MC, and χ V , whereas the second PCA component is mostly explained by the tortuosity and 3 Geofluids the third by the pore size and the porosity.95% of the variance can be explained by these first three PCA components.It is worth noting that, in contrast to the other geometric parameters, the micro-CT sandpack tortuosity is systematically lower than that of the synthetic media (Figures 1(k)-1(o)).Such a feature is not explained by the geometric characteristics captured by the variogram (or related to it), or by the first three Minkowski functionals.It is probably due to the limitations of variogram-based generators.As the name implies, such methods can generate random fields with a specified mean and variogram function, which are their one-point and two-point statistics, respectively.Higher order statistics can, by definition, not be reproduced, limiting the range of possible structures that such methods can produce (see, e.g., [47]).One particular feature that cannot be directly reproduced, and is often discussed in the literature, is the presence of long-ranging highlyconnected flow paths [48].Within the scope of our study, we do not consider this a major issue, since such longrange structures are hardly present at the pore scale.This   4 Geofluids notion is confirmed by our analysis, since the porous media generated using the Gaussian random fields are described effectively by the Minkovsky functionals and show properties comparable to the real media.

Lattice Boltzmann Method.
The lattice Boltzmann method is a computational fluid dynamics technique that solves fluid flow and transport processes within complex structures such as porous media [37,49,50].The mainstay of the LBM is the discretized Boltzmann equation (or lattice Boltzmann equation), which is based on kinetic theory and describes the streaming and collision of particles.The macroscopic variables are defined along each node of the computational domain as moments of the particle distribution  For diffusion simulations, no pressure gradient is applied and we assume for each run an initially homogenous concentration in a subdomain (i.e., C ini 50 ≤ x ≤ 250, 50 ≤ y ≤ 250, 50 ≤ z ≤ 250 = 1) and a concentration equals to 0 otherwise.A fixed null concentration is assigned at the inlet C inlet 0 ≤ x ≤ 300,0 ≤ y ≤ 300, z = 0 = 0 and the outlet C outlet 0 ≤ x ≤ 300,0 ≤ y ≤ 300, z = L = 0 of the domain, while periodic boundaries are set at the other faces (Figure 2(b)).
Finally, to simulate the combined solute advection and diffusion process, the solute transport simulation is initiated after the fluid flow reaches a steady state.The advection-diffusion simulations use the same initial conditions as the diffusion simulations, but a Dirichlet boundary condition for the concentration is assigned at the inlet face of the domain C inlet 0 ≤ x ≤ 300,0 ≤ y ≤ 300, z = 0 = 0 and a Neumann boundary condition ∂C/∂z = 0 with a secondorder finite difference scheme [52] is assigned to the outlet 0 ≤ x ≤ 300,0 ≤ y ≤ 300, z = L .Periodic boundary conditions are set as in the previous diffusive simulations (Figure 2(c)).
where ν is the kinematic viscosity of the fluid (assumed constant).Following Mostaghimi et al. (2016, [16]), we use the ratio π/SSA as the characteristic length scale.This choice is motivated since the specific surface area is readily computed for each medium.It represents a proper descriptor of the sample as it is directly related to the viscous forces acting at the fluid-rock interface.U z is the average seepage velocity along the main flow direction, calculated using the passing streamlines through the porous medium having a length that exceeds the domain side: where L s is the length of each streamline and N the number of streamlines passing through the porous media.The calculation of the seepage velocity suffers from the underestimation that may occur when the porosity is used in lieu of the effective porosity, because the latter is unknown [23].However, given the porosity range (between 0.3 and 0.4), the homogeneity, and relatively narrow pore size distribution of the media, we expect the effective and true porosity to closely match for each medium.Hereafter, the porosity will therefore be used in lieu of the effective porosity.
Finally, we use the original formula of Kozeny-Carman that provides a simple model to relate permeability to most of the descriptors discussed above: where β is a geometric factor (also called Kozeny constant) and depends on the shape of the cross sections of the channels.The original value of the Kozeny constant equals to 5, even if generally accepted literature values of β range from 2 for circular pores and 3 for flat cracks [46] to 10 [53].

Geometrical Control on Solute Diffusive Transport.
In the diffusion simulations, we compute for each time t the 5 Geofluids normalized solute mass within the computational spatial domain as follows: where This ratio quantifies the fraction of solute lost over time for the volume of porous medium initially occupied by the solute.By virtue of the initial conditions, this ratio is equal to 1 at t = 0 and decreases accordingly to the efficiency of the diffusive solute transport of each medium.In order to analyze the diffusion simulations for all porous media, we compared the time t * necessary for each sample to reach a selected cutoff of the initial mass, set at 10%.
The time t * is compared with the effective diffusion coefficient D eff , which is different from the molecular bulk diffusion coefficient.Several theoretical and empirical laws are proposed in the literature, mostly based on relationships that include either porosity and tortuosity combined or the tortuosity only [54][55][56][57].We test the following relationships: These different relationships vary depending on the spatial scale of investigation.For example, considering a macroscopic description of the effective diffusion coefficient requires both tortuosity and porosity.However, at the microscopic scale, the porosity is generally not considered [58].There is no clear agreement on linear or quadratic dependence on the tortuosity [55,56].The pressure gradient for all advection-diffusion simulations is fixed to 1.5 Pa/mm.Consequently, the Re number for each porous medium differs slightly because of the actual value of SSA and seepage velocity U z .The advection-diffusion simulations are carried out over a range of Pe values from 10 1 to 10 −2 for 9 samples selected over the initial 52: the two scanned samples and 7 synthetic media.
We calculate the hydrodynamic dispersion coefficient using the method of temporal moments [59][60][61][62].The dispersion coefficient is estimated from the first three temporal moments of the concentration breakthrough curves at the outlet of the computational domain.The n th temporal moment of a concentration distribution at the location z = L is defined as follows: We also define the second central moment as M 2c = M 2 − M 1 2 /M 0 .The longitudinal dispersion coefficient is given by 6 Geofluids Then, we compare the values of D L /D m computed through equation (11) with typical relationship used to estimate D L , as the following function of the Peclet number: where the exponent α varies from 1 to 2. Here, we used α = 1 2 accordingly to Mostaghimi et al. [16] and Icardi et al. [36].

Results and Discussion
In this section, the influence of the porous medium structure on fluid flow (Section 3.1) and solute diffusion (Section 3.2) is analyzed for all generated porous media as well as the two scanned sandpacks, while advection-diffusion simulations (Section 3.3) are described for nine samples selected among all porous media.
3.1.Fluid Flow Simulations.Figure 3 illustrates the pore-scale fluid flow for two of the generated porous media (Figures 3(a) and 3(b)) and the two sandpack samples (Figures 3(c) and 3(d)).The porous media EXX exhibit a very rough pore surface, whereas the GXX media show a smoother pore surface.The low velocities clearly define a creeping flow regime (Re <<1).For a detailed discussion about semivariogram functions adopted in the geostatistical model, we address the interested readers to Di Palma et al. [23].
Figure 4 shows the behavior of the seepage velocity with respect to the various geometric parameters analyzed.The seepage velocity relation to the pore size and porosity is approximately linear for each set of the two semivariogram models (Figures 4(a) and 4(b)).As pointed out in a previous study [23], the difference between the two semivariogram models is clear (velocities differ by two orders of magnitude), whereas the variations of porosity result in a vertical shift of the linear dispersion.The seepage velocity versus SSA displays a slightly nonlinear, inverse relation (Figure 4(c)), almost independent of the porosity.A 2nd-degree polynomial function has been fitted to highlight the relation between the seepage velocity and the SSA (the corresponding coefficient of determination R 2 is shown on the plot).The dependence of the seepage velocity on the MC (Figure 4(d)) appears clearer than its dependence on the SSA, the coefficient of determination being higher.This finding is quite intuitive, meaning how high velocities, which are responsible for the average flow rate, are controlled by the shape/size of channels passing through [26,28].
The trend between the Euler number and the seepage velocity (Figure 4(e)) shows a behavior similar to U z MC (Figure 4(d)).It means that the mean curvature and the Euler number provide similar information when describing the seepage velocity, even if the former is preferable because of its physical significance compared to the connectivity information provided by χ V that is purely topological.
Finally, the seepage velocity as a function of either tortuosity (Figure 4(f)) or porosity and tortuosity (Figures 4(g) and 4(h)) well identifies the generated porous media and the scanned samples in two subsets, although the values appear slightly scattered.
It is interesting to relate the computed seepage velocity to the permeability calculated using the Kozeny-Carman equation (equation ( 7)).The clear linear correlation among the generated media (i.e., black triangles and blue circles in Figure 5) suggests the proper combination of parameters governing the seepage velocity in the tested porous media.The wide range of permeability values also highlights the wide range explored through the two end-member semivariogram models.The values obtained for the sandpack samples (i.e., red squares) fall within the range of values obtained for the GXX porous media.The specific surface area represents the most sensitive parameter in the KC equation.Indeed, the porosities are identical for both semivariograms models; tortuosity varies over a limited range of values, whereas the SSA of EXX media is approximately two or three times greater than that of the GXX.
Because of the linearity between calculated seepage velocity and KC permeability, we can estimate the geometric factor (β) in equation ( 7) simply as the ratio of the imposed pressure gradient for all the fluid flow simulations over the slope of the fitting line.We obtain a Kozeny constant equals to 5.13, very similar to Carman's finding [46,63].
The analysis of the results for the advective flow shows that the SSA and the MC can be considered the most important/leading proxies of the advective flow.However, according to equation (7), the only specific surface area is not able to describe the linear relation between the seepage velocity and imposed pressure gradient (essentially the permeability), but it needs to be combined with the porosity and tortuosity to provide the best relationship (R 2 = 0 99).On the contrary, the permeability at the REV scale can be represented by the solely mean curvature, although with a small degradation of the performance (R 2 > 0 90).

Diffusion Simulations.
We analyze the behavior of the solute diffusive transport over time (equation ( 8)) and with the time t * as a metric for the efficiency of the diffusion process (Section 2.4.2).In analogy to the analysis of the seepage velocity, each Minkowski functional is plotted against t * (Figure 6).Figures 6(a)-6(e) did not show a clear correlation with t * .However, it is worth to note how the porosity is a good descriptor within each model adopted to generate porous media.
A common relationship can be observed when considering the inverse tortuosity (Figure 6(f)), with the exception of the porous media having the lowest correlation length (i.e., the upper right black triangles).This behavior results from a limited spatial resolution for the lowest correlation length.Excluding such realizations, a 2nd-degree polynomial function has been fitted and the corresponding coefficient of determination R 2 is shown on the plot to highlight the relation between the time t * and the inverse of the tortuosity squared.Finally, introducing the porosity in the descriptor (Figures 6(g) and 6(h)) splits the diffusive response accordingly to the two correlation models, in line with the results presented in Figure 6(a).In the following, the 1/τ 2 is considered as a unique proxy of the diffusive flow.

Advection-Diffusion
Simulations.The analysis of seepage velocity and time t * suggests that two pore-scale geometrical properties can be selected as proxies to describe either the advective or the diffusive transport (i.e., MC and 1/τ 2 ).This allows reducing the number of descriptors to only two.This two-parameter space is shown in Figure 7.It is worth noting that the selected media do not correspond exactly to the outermost points of this parameter space because the extreme values of generated porous media highlight different uncertainty issues: one related to low space resolution (i.e., EXX with lowest ξ corresponding to the highest MC values) and the second one to the unrepresentativeness of the domain length compared to the correlation length (i.e., GXX with highest ξ corresponding to the lowest MC values).
The results of the advection-diffusion simulations provide information that can be used to infer transport parameters at the REV scale.We estimate the hydrodynamic dispersion coefficient from the concentration breakthrough curves of each run using temporal moments (equation ( 11)) and compare them with published relationships (equation ( 12)).The Peclet number is estimated directly using the geometric properties of the porous media, where the seepage velocity is inferred from the KC formula, as validated in the previous section: Figure 8 shows the calculated longitudinal dispersion normalized by the molecular diffusivity for each simulation.We find a good match between the values calculated by the temporal moments and the relationship based on the geometric properties entering the definition of the Peclet number.Associated uncertainty has been estimated through the   [29,30,34].Moreover, the non-Fickian behavior is more evident under high advection-dominant transport [32].
One possible limit to the application of our findings to natural media lies in the assumption of geostatistical homogeneity of the generated porous media, which could prevent the extension of our results to porous media that exhibit multiple length scales and/or anisotropy.Future works should focus on the ability to generate realistic media displaying either or both heterogeneity and anisotropy.We do anticipate two different approaches to address this challenge.First, one could keep the present workflow but amend the procedure adopted in this paper by one additional transformation that allows for varying degrees of tortuosity (similar to the approach by [48]).Second, one could switch to the use of generators able to directly reproduce multipoint statistics [64,65].

Conclusions
This study investigates the impact of various porous medium geometric descriptors on the transport of non-reactive solutes.To that end, we performed pore-scale flow and transport simulations on 3D synthetic porous media and micro-CT images of sandpacks using the lattice Boltzmann method.Synthetic porous media are generated with a geostatistically  based approach, where the pore-scale characteristics are described by semivariogram models with a specified correlation length.We carried out a sensitivity analysis over different synthetic porous medium features to identify proxies that best describe advection and diffusive transport.Our findings suggest that the advective flow is primarily affected by the second and third Minkowski functionals (specific surface area and mean curvature).However, according to the Kozeny-Carman equation, macroscopic permeability necessitates three pore-scale parameters to be estimated (porosity, tortuosity, and specific surface area), while similar results can be obtained using the only mean curvature as proxy.
Results of the diffusion simulations show that the effective diffusion coefficient scales as the inverse of the tortuosity squared.
Finally, we tested the possibility to estimate the hydrodynamic dispersion coefficients using only the geometric properties of porous media and the applied pressure gradient.Such estimates from pore-scale geometric features have been compared with the coefficients calculated through the temporal moment method from the advection-diffusion simulations with Peclet numbers ranging from 10 -2 to 10 1 .We concluded, for the range of tested porous media, that an estimation of the hydrodynamic dispersion can be obtained from pore-scale properties: considering all tested Peclet numbers, the uncertainty of hydrodynamic dispersion is equal to 43%, when compared with that obtained from advection-diffusion simulations.Limiting such estimation to simulations characterized by Pe > 10 −1 , uncertainty decreases to 30%.11) and the dispersion coefficients obtained from Mostaghimi et al. [16] and Icardi et al. [36] and defined in equation ( 12).

Figure 1 :
Figure1: Relationships between values of correlation length (ξ), Minkowski functionals-porosity (φ), specific surface area (SSA), mean curvature (MC) and Euler number (χ V )-and tortuosity τ for generated (black triangles and blue circles) and scanned porous media (SP) as red squares used for numerical simulations.Porous media generated by exponential and Gaussian functions are named EXX and GXX, respectively, where XX refers to the correlation length of each medium.

2. 3 .
Numerical Simulation Setup.Numerical simulations for both flow and transport are carried out using Palabos[51], an open-source software for computational fluid dynamics based on the LBM.All calculations are performed on the cluster facility Euler at ETH Zurich.Initial and boundary conditions for all simulations are summarized in Figure2.Incompressible viscous flow is simulated directly through the pore space of the porous media by solving the Lattice Boltzmann equation for fluid flow.The fluid flow is driven by a fixed pressure gradient imposed along the z direction at the inlet (higher pressure at z = 0) and the outlet of the domain (lower pressure at z = L).Periodic boundary conditions are applied to the other edges of the computational domain.The solid-liquid interfaces along the grain boundaries are treated as a no-slip condition u = 0 , where u is the local fluid velocity.Simulations run until the fluid flow reaches a steady state (Figure2(a)).

2. 4 .
Quantitative Descriptors of the Simulated Processes 2.4.1.Descriptors of Flow Simulations.To describe the flow features, the Reynolds number is calculated as follows:

2. 4 . 3 .
Geometrical Control on Combined Advective-Diffusive Solute Transport.We perform the advection-diffusion simulations over a range of transport regimes, from diffusion to advection dominated.These simulations are characterized by the Peclet (Pe) number, comparing diffusive and advective flux and being defined as the product of Reynolds and Schmidt numbers Sc = υ/D m .We test three different Sc number values, i.e., 320, 32, and 3.2, obtained by varying the magnitude of the molecular diffusivity D m to avoid flow field variations due to kinematic viscosity of the fluid and to obtain faster simulations without affecting the transport process.

Figure 2 :
Figure 2: Schematic representation of the boundary conditions and initial conditions imposed for fluid flow (a), diffusion simulations (b), and advection-diffusion (c).

Figure 7 :
Figure 7: Parameter space of MC and 1/τ 2 for all generated porous media and CT image sandpack samples, obtained after fluid flow and diffusion simulations.Crossed markers indicate the sample selected for the advection-diffusion simulations.

Figure 8 :
Figure 8: Comparison between the hydrodynamic dispersion coefficients calculated from the advection-diffusion simulation using equation (11) and the dispersion coefficients obtained from Mostaghimi et al.[16] and Icardi et al.[36] and defined in equation(12).

Table 1 :
Characteristics of generated and scanned porous media used for numerical simulations.
which interact locally following simple microscopic laws.The LBM allows us to model several physical phenomena by means of different particle distribution functions, which are properly coupled to reproduce the process of interest (i.e., advection and diffusion within the single-phase fluid flow).Here, two particle distribution functions for fluid hydrodynamic f i x, t and concentration field g i x, t are used and a one-way coupling (the velocity of the fluid is used to calculate the advective term of equation S2.4 in Supplementary Materials) is adopted, since the transport of a non-reactive solute does not affect the fluid flow (Parmigiani et al., 2009).In LBM, all dimensional units are expressed in lattice units and, conventionally, lattice spacing and time step are unitary.More details about the lattice Boltzmann algorithm are given in Supplementary Material S1.