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

Di Palma, Paolo Roberto; Guyennon, Nicolas; Parmigiani, Andrea; Huber, Christian; Heβe, Falk; Romano, Emanuele

doi:https://doi.org/10.1155/2019/6810467

Geofluids

On this page

Abstract Introduction Materials and Methods Results and Discussion Conclusions Data Availability Additional Points Conflicts of Interest Supplementary Materials References Copyright Related Articles

Special Issue

Contribution of Pore-Scale Approach to Macroscale Geofluids Modelling in Porous Media

View this Special Issue

Research Article | Open Access

Volume 2019 | Article ID 6810467 | https://doi.org/10.1155/2019/6810467

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

Paolo Roberto Di Palma,¹Nicolas Guyennon,¹Andrea Parmigiani,²Christian Huber,³Falk Heβe,⁴and Emanuele Romano¹

Academic Editor: Andri Stefansson

Received17 Jul 2018

Revised30 Nov 2018

Accepted22 Jan 2019

Published03 Apr 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 10¹ to 10^-2.

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 (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 ()), 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 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 and , 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³, resulting in a number of voxels equal to 300³. 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 300³ 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 over the total volume of the medium ().

The second functional is the specific surface area, as the interstitial surface area of the voids and solids: where 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, 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 are the isolated objects, are interconnected regions or loops, and are the completely enclosed cavities. Broadly speaking, a negative corresponds to high pore space connectivity [45].

We calculate the tortuosity following the definition of hydraulic tortuosity in Clennell [46] as , being the average length of the fluid paths and is the straight-line length in the direction of flow [46]. 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 for a given set of starting points , where , , and 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 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 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 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 , 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 , 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 and concentration field 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 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 Figure 2.

(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 the direction at the inlet (higher pressure at ) and the outlet of the domain (lower pressure at ). 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 , where 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., ) and a concentration equals to 0 otherwise. A fixed null concentration is assigned at the inlet and the outlet 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 and a Neumann boundary condition with a second-order finite difference scheme [52] is assigned to the outlet . 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: where is the kinematic viscosity of the fluid (assumed constant). Following Mostaghimi et al. (2016, [16]), we use the ratio 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. 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 is the length of each streamline and 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].

2.4.2. Geometrical Control on Solute Diffusive Transport

In the diffusion simulations, we compute for each time the 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 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 necessary for each sample to reach a selected cutoff of the initial mass, set at 10%.

The time is compared with the effective diffusion coefficient , 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:

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 () number, comparing diffusive and advective flux and being defined as the product of Reynolds and Schmidt numbers . We test three different number values, i.e., 320, 32, and 3.2, obtained by varying the magnitude of the molecular diffusivity 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, the number for each porous medium differs slightly because of the actual value of SSA and seepage velocity . The advection-diffusion simulations are carried out over a range of values from 10¹ to 10⁻² 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 ^th temporal moment of a concentration distribution at the location is defined as follows:

We also define the second central moment as . The longitudinal dispersion coefficient is given by

Then, we compare the values of computed through equation (11) with typical relationship used to estimate , as the following function of the Peclet number: where the exponent varies from 1 to 2. Here, we used 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 (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 (). For a detailed discussion about semivariogram functions adopted in the geostatistical model, we address the interested readers to Di Palma et al. [23].

(a)

(b)

(c)

(d)

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 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 ()—and literature relationships for effective diffusivity, , , and , 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 () 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 (Figure 4(e)) shows a behavior similar to (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 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 (). 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 ().

3.2. Diffusion Simulations

We analyze the behavior of the solute diffusive transport over time (equation (8)) and with the time 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 (Figure 6). Figures 6(a)–6(e) did not show a clear correlation with . However, it is worth to note how the porosity is a good descriptor within each model adopted to generate porous media.

Figure 6

as a function of correlation length (), Minkowski functionals—porosity (φ), specific surface area (SSA), mean curvature (MC), and Euler number ()—and literature relationships for effective diffusivity, , , and , 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 () 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 (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 is shown on the plot to highlight the relation between the time 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 is considered as a unique proxy of the diffusive flow.

3.3. Advection-Diffusion Simulations

The analysis of seepage velocity and time 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 ). 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 symmetric mean absolute percentage error (SMAPE) defined as follows: resulting in a value equal to 43%. The uncertainty decreases for increasing (30% for ). 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 , 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.

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 10¹. 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 , uncertainty decreases to 30%.

Data Availability

The micro-CT images of sandpacks used to support the findings of this study have been deposited in the Figshare digital repository (https://www.imperial.ac.uk/earth-science/research/research-groups/perm/research/pore-scale-modelling/micro-ct-images-and-networks/). The generated porous media and simulation results used to generate figures are available upon request via e-mail to the first author.

Additional Points

Highlights. (i) We perform pore-scale flow simulations on synthetic porous media and micro-CT images. (ii) Advective transport is primarily affected by the second and third Minkowski functionals. (iii) Diffusive transport is a function of the inverse of the squared tortuosity. (iv) Hydrodynamic dispersion coefficients are estimated using geometric properties.

Conflicts of Interest

The authors declare that there is no conflict of interest regarding the publication of this paper.

Supplementary Materials

The details/description of the lattice Boltzmann method for single-phase fluid flow and non-reactive advection-diffusion simulations and the set of parameters used in this work. (Supplementary Materials)

References

J. W. Delleur, The Handbook of Groundwater Engineering, CRC press, 2006.
M. Rolle and P. K. Kitanidis, “Effects of compound-specific dilution on transient transport and solute breakthrough: a pore-scale analysis,” Advances in Water Resources, vol. 71, pp. 186–199, 2014.
View at: Publisher Site | Google Scholar
A. Guadagnini, M. J. Blunt, M. Riva, and B. Bijeljic, “Statistical scaling of geometric characteristics in millimeter scale natural porous media,” Transport in Porous Media, vol. 101, no. 3, pp. 465–475, 2014.
View at: Publisher Site | Google Scholar
Z. Liu, A. Herring, C. Arns, S. Berg, and R. T. Armstrong, “Pore-scale characterization of two-phase flow using integral geometry,” Transport in Porous Media, vol. 118, no. 1, pp. 99–117, 2017.
View at: Publisher Site | Google Scholar
C. Noiriel, P. Gouze, and B. Madé, “3D analysis of geometry and flow changes in a limestone fracture during dissolution,” Journal of Hydrology, vol. 486, pp. 211–223, 2013.
View at: Publisher Site | Google Scholar
T. D. Scheibe, W. A. Perkins, M. C. Richmond et al., “Pore-scale and multiscale numerical simulation of flow and transport in a laboratory-scale column,” Water Resources Research, vol. 51, no. 2, pp. 1023–1035, 2015.
View at: Publisher Site | Google Scholar
H. Andrä, N. Combaret, J. Dvorkin et al., “Digital rock physics benchmarks—part II: computing effective properties,” Computers & Geosciences, vol. 50, pp. 33–43, 2013.
View at: Publisher Site | Google Scholar
M. J. Blunt, B. Bijeljic, H. Dong et al., “Pore-scale imaging and modelling,” Advances in Water Resources, vol. 51, pp. 197–216, 2013.
View at: Publisher Site | Google Scholar
T. Bultreys, W. De Boever, and V. Cnudde, “Imaging and image-based fluid transport modeling at the pore scale in geological materials: a practical introduction to the current state-of-the-art,” Earth-Science Reviews, vol. 155, pp. 93–128, 2016.
View at: Publisher Site | Google Scholar
C. J. Werth, C. Zhang, M. L. Brusseau, M. Oostrom, and T. Baumann, “A review of non-invasive imaging methods and applications in contaminant hydrogeology research,” Journal of Contaminant Hydrology, vol. 113, no. 1-4, pp. 1–24, 2010.
View at: Publisher Site | Google Scholar
D. Wildenschild and A. P. Sheppard, “X-ray imaging and analysis techniques for quantifying pore-scale structure and processes in subsurface porous medium systems,” Advances in Water Resources, vol. 51, pp. 217–246, 2013.
View at: Publisher Site | Google Scholar
J. D. Hyman and C. L. Winter, “Stochastic generation of explicit pore structures by thresholding Gaussian random fields,” Journal of Computational Physics, vol. 277, pp. 16–31, 2014.
View at: Publisher Site | Google Scholar
J. T. Fredrich, A. A. DiGiovanni, and D. R. Noble, “Predicting macroscopic transport properties using microscopic image data,” Journal of Geophysical Research: Solid Earth, vol. 111, no. B3, 2006.
View at: Publisher Site | Google Scholar
M. Singh and K. K. Mohanty, “Permeability of spatially correlated porous media,” Chemical Engineering Science, vol. 55, no. 22, pp. 5393–5403, 2000.
View at: Publisher Site | Google Scholar
Y. Gao, X. Zhang, P. Rama et al., “Calculating the anisotropic permeability of porous media using the lattice Boltzmann method and X-ray computed tomography,” Transport in Porous Media, vol. 92, no. 2, pp. 457–472, 2012.
View at: Publisher Site | Google Scholar
P. Mostaghimi, B. Bijeljic, and M. Blunt, “Simulation of flow and dispersion on pore-space images,” SPE Journal, vol. 17, no. 4, pp. 1131–1141, 2012.
View at: Publisher Site | Google Scholar
P. Mostaghimi, M. J. Blunt, and B. Bijeljic, “Computations of absolute permeability on micro-CT images,” Mathematical Geosciences, vol. 45, no. 1, pp. 103–125, 2013.
View at: Publisher Site | Google Scholar
A. Ghassemi and A. Pak, “Pore scale study of permeability and tortuosity for flow through particulate media using Lattice Boltzmann method,” International Journal for Numerical and Analytical Methods in Geomechanics, vol. 35, no. 8, pp. 886–901, 2011.
View at: Publisher Site | Google Scholar
P. Lehmann, M. Berchtold, B. Ahrenholz et al., “Impact of geometrical properties on permeability and fluid phase distribution in porous media,” Advances in Water Resources, vol. 31, no. 9, pp. 1188–1204, 2008.
View at: Publisher Site | Google Scholar
A. Q. Raeini, M. J. Blunt, and B. Bijeljic, “Modelling two-phase flow in porous media at the pore scale using the volume-of-fluid method,” Journal of Computational Physics, vol. 231, no. 17, pp. 5653–5668, 2012.
View at: Publisher Site | Google Scholar
A. Q. Raeini, M. J. Blunt, and B. Bijeljic, “Direct simulations of two-phase flow on micro-CT images of porous media and upscaling of pore-scale forces,” Advances in Water Resources, vol. 74, pp. 116–126, 2014.
View at: Publisher Site | Google Scholar
D. Zhang, K. Papadikis, and S. Gu, “A lattice Boltzmann study on the impact of the geometrical properties of porous media on the steady state relative permeabilities on two-phase immiscible flows,” Advances in Water Resources, vol. 95, pp. 61–79, 2016.
View at: Publisher Site | Google Scholar
P. R. Di Palma, N. Guyennon, F. Heße, and E. Romano, “Porous media flux sensitivity to pore-scale geostatistics: a bottom-up approach,” Advances in Water Resources, vol. 102, pp. 99–110, 2017.
View at: Publisher Site | Google Scholar
P. K. Kang, P. de Anna, J. P. Nunes, B. Bijeljic, M. J. Blunt, and R. Juanes, “Pore-scale intermittent velocity structure underpinning anomalous transport through 3-D porous media,” Geophysical Research Letters, vol. 41, no. 17, pp. 6184–6190, 2014.
View at: Publisher Site | Google Scholar
M. Matyka, J. Gołembiewski, and Z. Koza, “Power-exponential velocity distributions in disordered porous media,” Physical Review E, vol. 93, no. 1, article 013110, 2016.
View at: Publisher Site | Google Scholar
M. Siena, M. Riva, J. D. Hyman, C. L. Winter, and A. Guadagnini, “Relationship between pore size and velocity probability distributions in stochastically generated porous media,” Physical Review E, vol. 89, no. 1, article 013018, 2014.
View at: Publisher Site | Google Scholar
M. Wang, Y. F. Chen, G. W. Ma, J. Q. Zhou, and C. B. Zhou, “Influence of surface roughness on nonlinear flow behaviors in 3D self-affine rough fractures: lattice Boltzmann simulations,” Advances in Water Resources, vol. 96, pp. 373–388, 2016.
View at: Publisher Site | Google Scholar
P. de Anna, B. Quaife, G. Biros, and R. Juanes, “Prediction of the low-velocity distribution from the pore structure in simple porous media,” Physical Review Fluids, vol. 2, no. 12, article 124103, 2017.
View at: Publisher Site | Google Scholar
B. Bijeljic, P. Mostaghimi, and M. J. Blunt, “Insights into non-Fickian solute transport in carbonates,” Water Resources Research, vol. 49, no. 5, pp. 2714–2728, 2013.
View at: Publisher Site | Google Scholar
B. Bijeljic, A. Raeini, P. Mostaghimi, and M. J. Blunt, “Predictions of non-Fickian solute transport in different classes of porous media using direct simulation on pore-scale images,” Physical Review E, vol. 87, no. 1, article 013011, 2013.
View at: Publisher Site | Google Scholar
M. Dentz, T. Le Borgne, A. Englert, and B. Bijeljic, “Mixing, spreading and reaction in heterogeneous media: a brief review,” Journal of Contaminant Hydrology, vol. 120-121, pp. 1–17, 2011.
View at: Publisher Site | Google Scholar
M. Dentz, M. Icardi, and J. J. Hidalgo, “Mechanisms of dispersion in a porous medium,” Journal of Fluid Mechanics, vol. 841, pp. 851–882, 2018.
View at: Publisher Site | Google Scholar
P. R. Di Palma, A. Parmigiani, C. Huber, N. Guyennon, and P. Viotti, “Pore-scale simulations of concentration tails in heterogeneous porous media,” Journal of Contaminant Hydrology, vol. 205, pp. 47–56, 2017.
View at: Publisher Site | Google Scholar
F. Gjetvaj, A. Russian, P. Gouze, and M. Dentz, “Dual control of flow field heterogeneity and immobile porosity on non-Fickian transport in Berea sandstone,” Water Resources Research, vol. 51, no. 10, pp. 8273–8293, 2015.
View at: Publisher Site | Google Scholar
D. L. Hochstetler, M. Rolle, G. Chiogna, C. M. Haberer, P. Grathwohl, and P. K. Kitanidis, “Effects of compound-specific transverse mixing on steady-state reactive plumes: insights from pore-scale simulations and Darcy-scale experiments,” Advances in Water Resources, vol. 54, pp. 1–10, 2013.
View at: Publisher Site | Google Scholar
M. Icardi, G. Boccardo, D. L. Marchisio, T. Tosco, and R. Sethi, “Pore-scale simulation of fluid flow and solute dispersion in three-dimensional porous media,” Physical Review E, vol. 90, no. 1, article 013032, 2014.
View at: Publisher Site | Google Scholar
S. Succi, The Lattice Boltzmann Equation: for Fluid Dynamics and Beyond, Oxford University Press, 2001.
H. J. Vogel, U. Weller, and S. Schlüter, “Quantification of soil structure based on Minkowski functions,” Computers & Geosciences, vol. 36, no. 10, pp. 1236–1245, 2010.
View at: Publisher Site | Google Scholar
F. Heße, V. Prykhodko, S. Schlüter, and S. Attinger, “Generating random fields with a truncated power-law variogram: a comparison of several numerical methods,” Environmental Modelling & Software, vol. 55, pp. 32–48, 2014.
View at: Publisher Site | Google Scholar
B. Matérn, Spatial Variation: Stochastic Models and Their Application to Some Problems in Forest Surveys and Other Sampling Investigations, Statens skogsforskningsinstitut, 1960.
B. Minasny and A. B. McBratney, “The Matérn function as a general model for soil variograms,” Geoderma, vol. 128, no. 3-4, pp. 192–207, 2005.
View at: Publisher Site | Google Scholar
O. Dubrule, “Indicator variogram models: do we have much choice?” Mathematical Geosciences, vol. 49, no. 4, pp. 441–465, 2017.
View at: Publisher Site | Google Scholar
F. D. E. Latief, B. Biswal, U. Fauzi, and R. Hilfer, “Continuum reconstruction of the pore scale microstructure for Fontainebleau sandstone,” Physica A: Statistical Mechanics and its Applications, vol. 389, no. 8, pp. 1607–1618, 2010.
View at: Publisher Site | Google Scholar
H. J. Vogel, “Topological characterization of porous media,” in Morphology of Condensed Matter, pp. 75–92, Springer, Berlin, Heidelberg, 2002.
View at: Publisher Site | Google Scholar
K. Wu, M. I. J. van Dijke, G. D. Couples et al., “3D stochastic modelling of heterogeneous porous media–applications to reservoir rocks,” Transport in Porous Media, vol. 65, no. 3, pp. 443–467, 2006.
View at: Publisher Site | Google Scholar
M. B. Clennell, “Tortuosity: a guide through the maze,” Geological Society, London, Special Publications, vol. 122, no. 1, pp. 299–344, 1997.
View at: Publisher Site | Google Scholar
N. Linde, P. Renard, T. Mukerji, and J. Caers, “Geological realism in hydrogeological and geophysical inverse modeling: a review,” Advances in Water Resources, vol. 86, pp. 86–101, 2015.
View at: Publisher Site | Google Scholar
B. Zinn and C. F. Harvey, “When good statistical models of aquifer heterogeneity go bad: a comparison of flow, dispersion, and mass transfer in connected and multivariate Gaussian hydraulic conductivity fields,” Water Resources Research, vol. 39, no. 3, 2003.
View at: Publisher Site | Google Scholar
T. Krüger, H. Kusumaatmaja, A. Kuzmin, O. Shardt, G. Silva, and E. M. Viggen, The Lattice Boltzmann Method: Principles and Practice, Springer, Cham, Switzerland, 2017.
View at: Publisher Site
M. C. Sukop and D. T. Thorne, Lattice Boltzmann Modeling: an Introduction for Geoscientists and Engineers, Springer, 2006.
View at: Publisher Site
J. Latt, “Palabos, parallel lattice Boltzmann solver,” 2009, January 2019, http://www.palabos.org.
View at: Google Scholar
C. Huber, B. Shafei, and A. Parmigiani, “A new pore-scale model for linear and non-linear heterogeneous dissolution and precipitation,” Geochimica et Cosmochimica Acta, vol. 124, pp. 109–130, 2014.
View at: Publisher Site | Google Scholar
A. Koponen, M. Kataja, and J. Timonen, “Permeability and effective porosity of porous media,” Physical Review E, vol. 56, no. 3, pp. 3319–3325, 1997.
View at: Publisher Site | Google Scholar
N. Epstein, “On tortuosity and the tortuosity factor in flow and diffusion through porous media,” Chemical Engineering Science, vol. 44, no. 3, pp. 777–779, 1989.
View at: Publisher Site | Google Scholar
B. Ghanbarian, A. G. Hunt, R. P. Ewing, and M. Sahimi, “Tortuosity in porous media: a critical review,” Soil Science Society of America Journal, vol. 77, no. 5, p. 1461, 2013.
View at: Publisher Site | Google Scholar
L. Shen and Z. Chen, “Critical review of the impact of tortuosity on diffusion,” Chemical Engineering Science, vol. 62, no. 14, pp. 3748–3755, 2007.
View at: Publisher Site | Google Scholar
J. Van Brakel and P. M. Heertjes, “Analysis of diffusion in macroporous media in terms of a porosity, a tortuosity and a constrictivity factor,” International Journal of Heat and Mass Transfer, vol. 17, no. 9, pp. 1093–1103, 1974.
View at: Publisher Site | Google Scholar
M. Matyka, A. Khalili, and Z. Koza, “Tortuosity-porosity relation in porous media flow,” Physical Review E, vol. 78, no. 2, article 026306, 2008.
View at: Publisher Site | Google Scholar
E. E. Adams and L. W. Gelhar, “Field study of dispersion in a heterogeneous aquifer: 2. Spatial moments analysis,” Water Resources Research, vol. 28, no. 12, pp. 3293–3307, 1992.
View at: Publisher Site | Google Scholar
O. A. Cirpka and P. K. Kitanidis, “Characterization of mixing and dilution in heterogeneous aquifers by means of local temporal moments,” Water Resources Research, vol. 36, no. 5, pp. 1221–1236, 2000.
View at: Publisher Site | Google Scholar
M. N. Goltz and P. V. Roberts, “Using the method of moments to analyze three-dimensional diffusion-limited solute transport from temporal and spatial perspectives,” Water Resources Research, vol. 23, no. 8, pp. 1575–1585, 1987.
View at: Publisher Site | Google Scholar
L. Pang, M. Goltz, and M. Close, “Application of the method of temporal moments to interpret solute transport with sorption and degradation,” Journal of Contaminant Hydrology, vol. 60, no. 1-2, pp. 123–134, 2003.
View at: Publisher Site | Google Scholar
P. C. Carman, “Fluid flow through granular beds,” Chemical Engineering Research and Design, vol. 75, pp. S32–S48, 1997.
View at: Publisher Site | Google Scholar
L. Y. Hu and T. Chugunova, “Multiple-point geostatistics for modeling subsurface heterogeneity: a comprehensive review,” Water Resources Research, vol. 44, no. 11, 2008.
View at: Publisher Site | Google Scholar
G. Mariethoz and S. Lefebvre, “Bridges between multiple-point geostatistics and texture synthesis: review and guidelines for future research,” Computers & Geosciences, vol. 66, pp. 66–80, 2014.
View at: Publisher Site | Google Scholar

Copyright

Copyright © 2019 Paolo Roberto Di Palma et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

PDF Download Citation

Download other formats

Order printed copies

Views

1534

Downloads

1625

Citations