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

Authors: Paolo Roberto Di Palma; Nicolas Guyennon; Andrea Parmigiani; Christian Huber; Falk Heβe; Emanuele Romano
Journal: Geofluids
 (2019)

---

## Abstract

Transport processes in porous media have been traditionally studied through the parameterization of macroscale properties, by means of volume-averaging or upscaling methods over a representative elementary volume. The possibility of upscaling results from pore-scale simulations, to obtain volume-averaging properties useful for practical purpose, can enhance the understanding of transport effects that manifest at larger scales. Several studies have been carried out to investigate the impact of the geometric properties of porous media on transport processes for solute species. However, the range of pore-scale geometric properties, which can be investigated, is usually limited to the number of samples acquired from microcomputed tomography images of real porous media. The present study takes advantage of synthetic porous medium generation to propose a systematic analysis of the relationships between geometric features of the porous media and transport processes through direct simulations of fluid flow and advection-diffusion of a non-reactive solute. Numerical simulations are performed with the lattice Boltzmann method on synthetic media generated with a geostatistically based approach. Our findings suggest that the advective transport is primarily affected by the specific surface area and the mean curvature of the porous medium, while the effective diffusion coefficient scales as the inverse of the tortuosity squared. Finally, the possibility of estimating the hydrodynamic dispersion coefficient knowing only the geometric properties of porous media and the applied pressure gradient has been tested, within the range of tested porous media, against advection-diffusion simulations at low Reynolds (<10-1) and Peclet numbers ranging from 101 to 10-2.


---

## Body

## 1. 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–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–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 well-established use of pore-scale approach consists of estimating averaged properties, such as permeability through direct single-phase fluid flow simulations [6, 13–18] and relative permeability for multiphase fluid flows [4, 19–22]. Recently, some authors focused on the relationship between velocity distribution and pore structure [23–28] as well as transport processes [2, 24, 29–36].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].


## 2. Materials and Methods

To identify the dominant pore-scale parameters that control the transport of non-reactive solutes, LBM simulations (Section2.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 MaterialS3.


### 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 mm3, resulting in a number of voxels equal to 3003. 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 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 3003 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:(1)ϕ=VvoidsV,defined as the volume of the voids Vvoids 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:(2)SSA=AsV,where As 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:(3)MC=W1V.Here,W1 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:
(4)χV=N−C+HV,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 τ=Le/L, Le being the average length of the fluid paths and L is the straight-line length in the direction of flow [46]. Le 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 xṫ,yṫ,zṫ=uxt,yt,zt,vxt,yt,zt,wxt,yt,zt for a given set of starting points xt=0,yt=0,zt=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 Table1, together with the corresponding estimated geometric properties for the scanned samples.Table 1
Characteristics of generated and scanned porous media used for numerical simulations.
Exponential model media
Gaussian model media
X-ray scanned sandpack media
ID
E09
E15
E21
E27
E33
G09
G15
G21
G27
G33
F42C
LV60A
ξ (mm)
0.1
0.2
0.2
0.3
0.3
0.1
0.2
0.2
0.3
0.3
0.16
0.125
φ (-)
0.30/0.325/0.35/0.375/0.40
0.30/0.325/0.35/0.375/0.40
0.33
0.38Figure1 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.Figure 1
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.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 Figure1 (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 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 highly-connected flow paths [48]. Within the scope of our study, we do not consider this a major issue, since such long-range structures are hardly present at the pore scale. This 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.


### 2.2. 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 function fix,t, 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 fix,t and concentration field gix,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.


### 2.3. Numerical Simulation Setup

Numerical simulations for both flow and transport are carried out usingPalabos [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 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).
(a)
(b)
(c)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 thez 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 (Figure 2(a)).For diffusion simulations, no pressure gradient is applied and we assume for each run an initially homogenous concentration in a subdomain (i.e.,Cini50≤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 Cinlet0≤x≤300,0≤y≤300,z=0=0 and the outlet Coutlet0≤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 domainCinlet0≤x≤300,0≤y≤300,z=0=0 and a Neumann boundary condition ∂C/∂z=0 with a second-order 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)).


### 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:(5)Re=U¯z·π/SSAν,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:
(6)U¯z=∑s=1N∑i=1Lsuzxi,yi,ziN,where Ls 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:(7)k=ϕ3βτ2SSA2,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].


#### 2.4.2. Geometrical Control on Solute Diffusive Transport

In the diffusion simulations, we compute for each timet the normalized solute mass within the computational spatial domain as follows:
(8)MintM0=∑i=1Nx∑j=1Ny∑k=1NzCxi,yj,zk,tM0,where M0=∑i=1Nx∑j=1Ny∑k=1NzCxi,yj,zk,t=0. 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 timet∗ is compared with the effective diffusion coefficient Deff, 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–57]. We test the following relationships:
(9)Deff=1τ2Dm,ϕτDm,ϕτ2Dm.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].


#### 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=υ/Dm. We test three different Sc number values, i.e., 320, 32, and 3.2, obtained by varying the magnitude of the molecular diffusivity Dm to avoid flow field variations due to kinematic viscosity of the fluid and to obtain faster simulations without affecting the transport process.The pressure gradient for all advection-diffusion simulations is fixed to 1.5 Pa/mm. Consequently, theRe 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 101 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–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 nth temporal moment of a concentration distribution at the location z=L is defined as follows:
(10)Mn=∫0∞tncx,y,z=L,tdt.We also define the second central moment asM2c=M2−M12/M0. The longitudinal dispersion coefficient is given by
(11)DL=z2M2cM022M13.Then, we compare the values ofDL/Dm computed through equation (11) with typical relationship used to estimate DL, as the following function of the Peclet number:
(12)DLDm=1τ2+Peα,where the exponent α varies from 1 to 2. Here, we used α=1.2 accordingly to Mostaghimi et al. [16] and Icardi et al. [36].


## 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 mm3, resulting in a number of voxels equal to 3003. 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 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 3003 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:(1)ϕ=VvoidsV,defined as the volume of the voids Vvoids 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:(2)SSA=AsV,where As 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:(3)MC=W1V.Here,W1 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:
(4)χV=N−C+HV,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 τ=Le/L, Le being the average length of the fluid paths and L is the straight-line length in the direction of flow [46]. Le 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 xṫ,yṫ,zṫ=uxt,yt,zt,vxt,yt,zt,wxt,yt,zt for a given set of starting points xt=0,yt=0,zt=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 Table1, together with the corresponding estimated geometric properties for the scanned samples.Table 1
Characteristics of generated and scanned porous media used for numerical simulations.
Exponential model media
Gaussian model media
X-ray scanned sandpack media
ID
E09
E15
E21
E27
E33
G09
G15
G21
G27
G33
F42C
LV60A
ξ (mm)
0.1
0.2
0.2
0.3
0.3
0.1
0.2
0.2
0.3
0.3
0.16
0.125
φ (-)
0.30/0.325/0.35/0.375/0.40
0.30/0.325/0.35/0.375/0.40
0.33
0.38Figure1 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.Figure 1
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.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 Figure1 (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 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 highly-connected flow paths [48]. Within the scope of our study, we do not consider this a major issue, since such long-range structures are hardly present at the pore scale. This 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.


## 2.2. 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 function fix,t, 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 fix,t and concentration field gix,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.


## 2.3. Numerical Simulation Setup

Numerical simulations for both flow and transport are carried out usingPalabos [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 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).
(a)
(b)
(c)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 thez 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 (Figure 2(a)).For diffusion simulations, no pressure gradient is applied and we assume for each run an initially homogenous concentration in a subdomain (i.e.,Cini50≤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 Cinlet0≤x≤300,0≤y≤300,z=0=0 and the outlet Coutlet0≤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 domainCinlet0≤x≤300,0≤y≤300,z=0=0 and a Neumann boundary condition ∂C/∂z=0 with a second-order 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)).


## 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:(5)Re=U¯z·π/SSAν,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:
(6)U¯z=∑s=1N∑i=1Lsuzxi,yi,ziN,where Ls 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:(7)k=ϕ3βτ2SSA2,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].


### 2.4.2. Geometrical Control on Solute Diffusive Transport

In the diffusion simulations, we compute for each timet the normalized solute mass within the computational spatial domain as follows:
(8)MintM0=∑i=1Nx∑j=1Ny∑k=1NzCxi,yj,zk,tM0,where M0=∑i=1Nx∑j=1Ny∑k=1NzCxi,yj,zk,t=0. 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 timet∗ is compared with the effective diffusion coefficient Deff, 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–57]. We test the following relationships:
(9)Deff=1τ2Dm,ϕτDm,ϕτ2Dm.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].


### 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=υ/Dm. We test three different Sc number values, i.e., 320, 32, and 3.2, obtained by varying the magnitude of the molecular diffusivity Dm to avoid flow field variations due to kinematic viscosity of the fluid and to obtain faster simulations without affecting the transport process.The pressure gradient for all advection-diffusion simulations is fixed to 1.5 Pa/mm. Consequently, theRe 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 101 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–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 nth temporal moment of a concentration distribution at the location z=L is defined as follows:
(10)Mn=∫0∞tncx,y,z=L,tdt.We also define the second central moment asM2c=M2−M12/M0. The longitudinal dispersion coefficient is given by
(11)DL=z2M2cM022M13.Then, we compare the values ofDL/Dm computed through equation (11) with typical relationship used to estimate DL, as the following function of the Peclet number:
(12)DLDm=1τ2+Peα,where the exponent α varies from 1 to 2. Here, we used α=1.2 accordingly to Mostaghimi et al. [16] and Icardi et al. [36].


## 2.4.1. Descriptors of Flow Simulations

To describe the flow features, the Reynolds number is calculated as follows:(5)Re=U¯z·π/SSAν,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:
(6)U¯z=∑s=1N∑i=1Lsuzxi,yi,ziN,where Ls 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:(7)k=ϕ3βτ2SSA2,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].


## 2.4.2. Geometrical Control on Solute Diffusive Transport

In the diffusion simulations, we compute for each timet the normalized solute mass within the computational spatial domain as follows:
(8)MintM0=∑i=1Nx∑j=1Ny∑k=1NzCxi,yj,zk,tM0,where M0=∑i=1Nx∑j=1Ny∑k=1NzCxi,yj,zk,t=0. 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 timet∗ is compared with the effective diffusion coefficient Deff, 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–57]. We test the following relationships:
(9)Deff=1τ2Dm,ϕτDm,ϕτ2Dm.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].


## 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=υ/Dm. We test three different Sc number values, i.e., 320, 32, and 3.2, obtained by varying the magnitude of the molecular diffusivity Dm to avoid flow field variations due to kinematic viscosity of the fluid and to obtain faster simulations without affecting the transport process.The pressure gradient for all advection-diffusion simulations is fixed to 1.5 Pa/mm. Consequently, theRe 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 101 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–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 nth temporal moment of a concentration distribution at the location z=L is defined as follows:
(10)Mn=∫0∞tncx,y,z=L,tdt.We also define the second central moment asM2c=M2−M12/M0. The longitudinal dispersion coefficient is given by
(11)DL=z2M2cM022M13.Then, we compare the values ofDL/Dm computed through equation (11) with typical relationship used to estimate DL, as the following function of the Peclet number:
(12)DLDm=1τ2+Peα,where the exponent α varies from 1 to 2. Here, we used α=1.2 accordingly to Mostaghimi et al. [16] and Icardi et al. [36].


## 3. Results and Discussion

In this section, the influence of the porous medium structure on fluid flow (Section3.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

Figure3 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 3
Visualization of porous medium structures and fluid flow streamlines. (a) and (b) are two generated media withξ=0.21 mm and φ=0.35, whereas (c) and (d) are the scanned sandpacks with estimated ξ=0.16 mm and φ=0.33 (F42C) and ξ=0.125 mm and φ=0.38 (LV60A).
(a)
(b)
(c)
(d)Figure4 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 R2 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].Figure 4
Seepage velocity as a function of correlation length (ξ), Minkowski functionals—porosity (φ), specific surface area (SSA), mean curvature (MC), and Euler number (χV)—and literature relationships for effective diffusivity, 1/τ2, τ2/τ, and ϕ/τ2, respectively, for generated porous media (black triangles and blue circles for exponential (EXX) and Gaussian (GXX) variogram models, respectively) and scanned sandpacks (SP) (red squares). The coefficient of determination (R2) associated with the fitting functions for the SSA and MC (red lines) is reported in (c, d).The trend between the Euler number and the seepage velocity (Figure4(e)) shows a behavior similar to U¯zMC (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 (Figure4(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.Figure 5
Seepage velocity as function of permeability. Porous media generated by exponential and Gaussian variogram models are named EXX and GXX (black triangles and blue circles, respectively).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 (R2=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 (R2>0.90).


### 3.2. 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.Figure 6
t∗ as a function of correlation length (ξ), Minkowski functionals—porosity (φ), specific surface area (SSA), mean curvature (MC), and Euler number (χV)—and literature relationships for effective diffusivity, 1/τ2, ϕ/τ, and ϕ/τ2, respectively, for generated porous media (black triangles and blue circles for exponential (EXX) and Gaussian (GXX) variogram models, respectively) and scanned sandpacks (SP) (red squares). The coefficient of determination (R2) associated to the fitting functions for the MC (red lines) is reported in (f).A common relationship can be observed when considering the inverse tortuosity (Figure6(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 R2 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.


### 3.3. Advection-Diffusion Simulations

The analysis of seepage velocity and timet∗ 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).Figure 7
Parameter space of MC and1/τ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.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:
(13)Pe=U¯zπ/SSADm=πϕ2βτ2SSA3DmμdPL.Figure8 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 symmetric mean absolute percentage error (SMAPE) defined as follows:
(14)SMAPE=100n∑i=1nX−YX+Y/2,resulting in a value equal to 43%. The uncertainty decreases for increasing Pe (30% for Pe>10−1). In case of diffusion-dominant process, the mechanical dispersion vanishes and the effective diffusivity can be calculated by fitting the concentration breakthrough curves with analytical expressions. The shift from high to low Pe, at the same conditions of imposed pressure gradient, is a straightforward consequence of the SSA range of the porous media, highlighting again the impact of this geometric feature.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).The good estimation of the longitudinal dispersion through the temporal moments confirms the validity of Fickian transport regime for the investigated porous media, as also reported in Bijeljic et al. [30] and Icardi et al. [36] that used realistic packed geometries comparable to our synthetic media. The non-Fickian behavior shows up in more complex porous media, such as carbonate and sandstone rock, in which the solute experiences a very wide range of transit times across pores of a very different size [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].


## 3.1. Fluid Flow Simulations

Figure3 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 3
Visualization of porous medium structures and fluid flow streamlines. (a) and (b) are two generated media withξ=0.21 mm and φ=0.35, whereas (c) and (d) are the scanned sandpacks with estimated ξ=0.16 mm and φ=0.33 (F42C) and ξ=0.125 mm and φ=0.38 (LV60A).
(a)
(b)
(c)
(d)Figure4 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 R2 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].Figure 4
Seepage velocity as a function of correlation length (ξ), Minkowski functionals—porosity (φ), specific surface area (SSA), mean curvature (MC), and Euler number (χV)—and literature relationships for effective diffusivity, 1/τ2, τ2/τ, and ϕ/τ2, respectively, for generated porous media (black triangles and blue circles for exponential (EXX) and Gaussian (GXX) variogram models, respectively) and scanned sandpacks (SP) (red squares). The coefficient of determination (R2) associated with the fitting functions for the SSA and MC (red lines) is reported in (c, d).The trend between the Euler number and the seepage velocity (Figure4(e)) shows a behavior similar to U¯zMC (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 (Figure4(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.Figure 5
Seepage velocity as function of permeability. Porous media generated by exponential and Gaussian variogram models are named EXX and GXX (black triangles and blue circles, respectively).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 (R2=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 (R2>0.90).


## 3.2. 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.Figure 6
t∗ as a function of correlation length (ξ), Minkowski functionals—porosity (φ), specific surface area (SSA), mean curvature (MC), and Euler number (χV)—and literature relationships for effective diffusivity, 1/τ2, ϕ/τ, and ϕ/τ2, respectively, for generated porous media (black triangles and blue circles for exponential (EXX) and Gaussian (GXX) variogram models, respectively) and scanned sandpacks (SP) (red squares). The coefficient of determination (R2) associated to the fitting functions for the MC (red lines) is reported in (f).A common relationship can be observed when considering the inverse tortuosity (Figure6(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 R2 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.


## 3.3. Advection-Diffusion Simulations

The analysis of seepage velocity and timet∗ 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).Figure 7
Parameter space of MC and1/τ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.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:
(13)Pe=U¯zπ/SSADm=πϕ2βτ2SSA3DmμdPL.Figure8 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 symmetric mean absolute percentage error (SMAPE) defined as follows:
(14)SMAPE=100n∑i=1nX−YX+Y/2,resulting in a value equal to 43%. The uncertainty decreases for increasing Pe (30% for Pe>10−1). In case of diffusion-dominant process, the mechanical dispersion vanishes and the effective diffusivity can be calculated by fitting the concentration breakthrough curves with analytical expressions. The shift from high to low Pe, at the same conditions of imposed pressure gradient, is a straightforward consequence of the SSA range of the porous media, highlighting again the impact of this geometric feature.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).The good estimation of the longitudinal dispersion through the temporal moments confirms the validity of Fickian transport regime for the investigated porous media, as also reported in Bijeljic et al. [30] and Icardi et al. [36] that used realistic packed geometries comparable to our synthetic media. The non-Fickian behavior shows up in more complex porous media, such as carbonate and sandstone rock, in which the solute experiences a very wide range of transit times across pores of a very different size [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].


## 4. 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 101. 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%.



---

*Source: 6810467-2019-04-03.xml*