Abstract
We compute effective properties (i.e., permeability, hydraulic tortuosity, and diffusive tortuosity) of three different digital porous media samples, including in-line array of uniform shapes, staggered-array of squares, and randomly distributed squares. The permeability and hydraulic tortuosity are computed by solving a set of rescaled Stokes equations obtained by homogenization, and the diffusive tortuosity is computed by solving a homogenization problem given for the effective diffusion coefficient that is inversely related to diffusive tortuosity. We find that hydraulic and diffusive tortuosity can be quantitatively different by up to a factor of ten in the same pore geometry, which indicates that these tortuosity terms cannot be used interchangeably. We also find that when a pore geometry is characterized by an anisotropic permeability, the diffusive tortuosity (and correspondingly the effective diffusion coefficient) can also be anisotropic. This finding has important implications for reservoir-scale modeling of flow and transport, as it is more realistic to account for the anisotropy of both the permeability and the effective diffusion coefficient.
1. Introduction
When modeling subsurface flow and transport in a reservoir, it is common practice to represent the geological formation with effective properties instead of resolving the precise location of fluid and solid phase. Effective properties such as permeability, diffusivity, and tortuosity can be computed from pore-scale modeling of fluid flow and transport through a porous media sample (i.e., a core sample) taken from a geological formation. Permeability and effective diffusivity are the well-known and important properties used in Darcy- and reservoir-scale flow and transport equations, while tortuosity is given less importance and subject to multiple definitions. Tortuosity is generally defined as a ratio between the effective path traveled by a species and the unit length of the domain. As noted in [1], different types of tortuosity (geometric, hydraulic, diffusive, dispersive, and electrical) are distinguished by the type of species transport under consideration. Hydraulic tortuosity refers to the effective path traveled by a fluid particle driven by a force, while diffusive tortuosity refers to the effective path traveled by a species driven by molecular diffusion [1].
Pore-scale modeling is performed on digital samples of porous media, which can be obtained using imaging and conversion techniques as those described in [2, 3], or by using an algorithm that generates samples which possess the same statistical characteristics of a real porous media sample [4–6]. Table 1 presents a summary of studies which used either real or computer-generated porous media samples; real refers to a digital sample obtained by imaging and conversion to binary form, and computer-generated refers to a digital sample obtained by an algorithm or random reconstruction. We note that Table 1 is only a short list of recent works on the topic, and that pore-scale simulation to obtain effective properties dates back to Cancelliere et al. [7] and possibly earlier.
In regard to obtaining effective properties at the pore-scale, we are motivated to distinguish between two different types of tortuosity (namely, hydraulic and diffusive) that have been identified. While a few studies have compared tortuosity definitions and have shown they are indeed quantitatively different within the same pore structure [8, 9], other studies continue to use these definitions interchangeably or do not distinguish the type of tortuosity they are computing. The usefulness of hydraulic tortuosity becomes apparent when formulating a prediction for permeability, such as the Kozeny-Carman equation. Also, tortuosity can be seen as a more intuitive quantity for fluid flow and transport through a pore geometry than permeability or the effective diffusion coefficient.
The outline of this work is as follows: first we will conduct a literature review on the methods commonly used to compute the effective properties of digital porous media; then we will present the details of these methods we employed herein, followed by three example geometries (in-line array of uniform shapes, staggered-array of squares, and randomly distributed squares) and our obtained trends between porosity and the effective properties. We will conclude with the main findings of this work.
2. Literature Review on Effective Properties
2.1. Permeability
From the list given in Table 1, an important and commonly computed effective property is permeability. The notion of a porous media’s permeability can be understood in terms of Darcy’s equation, which computes a macroscale fluid velocity bywhere is the permeability tensor. The other quantities are fluid viscosity and density and , macroscale pressure , and gravitational acceleration . Eqn. (1) can be obtained mathematically by the method of homogenization [10, 11], which is a multiscale expansion of the fluid flow equations at the pore-scale (i.e., Stokes equations); see Section 3.
Out of the papers listed in Table 1, the general approach to compute the permeability field of the porous media sample is as follows:(1)Obtain pore-scale geometries: converted from real media sample images or generated by an algorithm that may or may not use media parameters such as grain size distribution and so on.(2)Solve Stokes flow in pore space of media: different solvers have been used (finite difference [12], finite element, lattice Boltzmann [13–16], explicit jump method [15], etc.), from academically developed codes, commercial software, and open-source software.(3)Compute permeability: assuming Darcy’s equation is valid for the pore structure and flow field results, permeability is computed based on volume averages of velocity and pressure field.
A further step is to fit the permeability data to a formula which proposes a relationship between permeability and other porous media properties, known as the Kozeny-Carman (KC) equation [17, 18].
2.2. Kozeny-Carman Equation
An idealized pore geometry of parallel cylindrical channels was used to formulate the Kozeny-Carman (KC) equation [17]. This equation proposes a relationship between porosity and permeability , as well as other parameters such as hydraulic tortuosity and specific surface area , and iswhere is a fitting parameter (sometimes referred to as the shape factor or the KC constant, though inconsistently) which is used to account for different channel configurations. A range of shape factors have been presented in literature [17, 18] which correspond to cylindrical, elliptical, rectangular, and triangular cross-sectional channel shapes. The parameter in (2) is the ratio between effective path traveled by the fluid particle and the length of the sample: . In order to be consistent with the literature we are referencing and comparing our results against, we call this the hydraulic tortuosity, although we recognize that other work refers to as the hydraulic tortuosity as similarly pointed out in [19]. The specific surface area is the ratio of fluid-solid interface area to the total volume in the domain:Notice that (2) is not written as direction-specific. In fact, Carrier III [20] mentions that a limitation of the KC equation is that it does not explicitly account for anisotropy (a direction-specific property) even though the permeability of a real geological formation is typically greater in the horizontal than in the vertical direction. Despite this limitation, the Kozeny-Carman equation could still apply to direction-specific flows, where and are the permeability and hydraulic tortuosity in the principal directions, respectively.
Various studies [14, 21, 22] include attempts at establishing a modification to the KC equation, or propose a shape factor for a specific class of porous media. Xu and Yu [23] compiled a list of 11 studies that proposed a modification to the KC equation, and each work corresponded to a different media type (textile, glass and fiber, square particles, sandstone, etc.). These proposed KC modifications included different parameters like effective porosity, percolation threshold, grain radius, fractal dimension, interconnectivity parameter, and others.
2.3. Tortuosity
2.3.1. Different Types of Tortuosity
As stated in the introduction, different types of tortuosities have been defined according to the type of species transport under consideration. Previous works have recognized this and some provided a comparison between specific tortuosity types [1, 8, 13]. What is important is that these forms of tortuosity should not be expected to be the same in the same porous structure [24], unless the pore size distribution is very narrow as pointed out in Ghanbarian et al. [8]. Perhaps unintentionally, a few previous studies have failed to distinguish the difference between tortuosity types. For example, Ohkubo [25] measured diffusive tortuosity through porous media that was represented by plate-like obstacles and used the diffusive tortuosity results to compute fitting coefficients present in Koponen et al.’s [26] tortuosity-porosity trend equation. However, Koponen et al.’s [26] trend was empirically derived from simulation results that computed hydraulic tortuosity. Also, Sun et al. [27] applied the method of homogenization to the diffusive transport equation in order to compute the effective diffusion coefficient of periodic unit cells and the corresponding tortuosity. What they did not point out was that their approach gave them the diffusive tortuosity, while some of the trends they compared their results to were with respect to hydraulic tortuosity. Unless tortuosity types are quantitatively identical in the same class of porous media, use of a trend that is specific to one type of tortuosity should not immediately be used to describe the trend of another type of tortuosity.
A few studies [8, 9, 28] have focused on quantifying the difference between two specific tortuosities in the same structure. For example, hydraulic and electrical tortuosity were compared in both David [28] and Zhang and Knackstedt [9]. David [28] used a network as an analog for the media pore space, and the study focused on quantifying the tortuosity of the critical or preferential pathway, not the tortuosity of the entire flow field. Zhang and Knackstedt [9] computed tortuosity using the entire flow field which they obtained by numerical simulations for fluid flow (via lattice Gas Automata) and electrical current (via finite difference). Both of these works found that the hydraulic tortuosity was higher than electrical tortuosity in the same structure (up to an order of magnitude in low porosity configurations [9]). Once again, it was emphasized in Ghanbarian et al.’s [8] review that due to the difference between tortuosity types their models are not interchangeable.
Regarding diffusive tortuosity, both Quintard et al. [29] and Valdés-Parada et al. [1] stated that its quantification becomes more complex when more than just passive diffusion is occurring in the system. Valdés-Parada et al. [1] showed how the quantification of diffusive tortuosity can be different depending on the consideration of different mass transport terms (i.e., passive diffusion only, or diffusion with convection, reaction, or hydrodynamic dispersion). They concluded that the consideration of passive diffusion leads to the only appropriate definition of tortuosity. Their work focused on mass transport of species through fluid and did not measure hydraulic tortuosity.
2.3.2. Various Tortuosity-Porosity Trends
Many authors have theoretically or empirically derived tortuosity as a function of porosity (see Tables 2 and 3, where the tortuosity corresponds to fluid flow (hydraulic) and diffusive transport, resp.). Similar tables that summarize the theoretically and empirically derived trends proposed in literature appear in Shen and Chen [30], Boudreau [31], and Ahmadi et al. [32]. The most recent and complete comparison of tortuosity trends (or models) that we have seen to date was made by Ghanbarian et al. [8], which includes an explicit comparison between many different tortuosity types, namely, geometric, hydraulic, electrical, diffusive, and even tortuosity models, for unsaturated porous media.
In the same way that the KC equation is specific to a given class of porous media, any porosity-tortuosity trend that is empirically derived depends on the porous media considered. Some authors rely on synthetic porous media, such as 2D randomly distributed squares [13, 19, 24, 26, 33, 34], 3D high porosity plate-like obstacles [25], or unit cells of centered uniform shapes, while others rely on real samples of media taken by imaging (via scanning electron microscopy (SEM) or X-ray microtomography) and digitization.
Additionally, it is important to note that empirically derived trends are applicable to a limited range of porosity only, as the trend was developed given simulated data of a certain porosity range. The porosity limit could be a result of the idealized geometry under consideration. Ghanbarian et al. [8] noticed that most tortuosity models were derived by focusing on the higher porosity structures and recommended that research moves more to focus on the lower porosity structure where the porosity approaches the percolation threshold (i.e., the minimum porosity for which flow is still possible through structure).
While many trends (or models) have been proposed to compute tortuosity as a function of porosity and other fitting parameters, Valdés-Parada et al. [1] stated that any definition of tortuosity should not be considered as a function of porosity but rather a function of the pore geometry only. The idea behind this statement seems to rest in the fact that there is no universally agreed upon trend between tortuosity and porosity as pointed out by Matyka et al. [24] and that the geometrical features of a pore structure should have more influence on the overall tortuosity of flow or transport through the geometry than the value of porosity. However, since it is possible that a unique relationship between tortuosity and porosity can indeed exist for special classes of porous media [24], the development of tortuosity-porosity trends remains to be an insightful topic of research.
Also, Liu and Kitanidis [35] stated that tortuosity is a tensorial property. Despite this statement, the majority of work does not report directional-specific quantities which may leave the reader with the impression that tortuosity is a scalar quantity. A scalar quantity could be appropriate in pore geometries that exhibit an anisotropic microstructure while still exhibiting an isotropic macrostructure. For example, pore geometry composed of randomly distributed squares is geometrically anisotropic at the microscale; however the effective properties computed for the REV can be isotropic in nature. In this work, we demonstrate this point by example.
3. Methods to Compute Effective Properties
3.1. Permeability
The method of homogenization has been used to mathematically derive Darcy’s equation from the Stokes equations. Through this derivation, the macroscopic property known as permeability is defined as the spatial average of a rescaled pore space velocity . In the following, we present the basic steps that lead to such a definition, and the reader is referred to [11] for more derivation details.
The problem is first defined by the macroscale (a) and microscale (b) systems shown in Figure 1, where the macroscale domain is comprised of repeating unit cells of size , which is comprised of both fluid and solid space. Due to periodicity, the macroscale domain of porosity can be represented by the unit cell with periodic boundaries. Assuming slow and steady flow of incompressible fluid through the pore space of the media, the governing equations are given by the Stokes equations as follows:where and and repeated indices imply summation, is pressure, is fluid viscosity, is the pore space fluid velocity, is fluid density, and is gravitational acceleration. The boundary condition at the fluid-solid interface is given bywhich comes from the no-flow and no-slip boundary conditions. Although the slip boundary condition is necessary in certain contexts (such as gas flow and non-Newtonian fluid), the no-slip boundary condition is a valid assumption for viscous flow in porous media and is widely used in modeling studies of pore-scale flow [12, 15, 24].
In the Stokes equations, and can vary within the pore space and thus are functions of the macroscale and microscale (); however , , and are treated as constants and are thus functions of the macroscale only. By defining a fast variable , the cell of unit length is scaled into the -cell (shown in Figure 1(b)), where . A multiscale expansion is used to write and in terms of their leading and higher order terms as follows:
Upon substitution of the expanded variables into (4)–(6) and then the collection of like order terms, the second-order pressure term and the third-order velocity term aresuch that and satisfywith the fluid-solid boundary condition given byIn the above formulation, and are functions of the microscale only and are independent of macroscale quantities , , and . They can be thought of as rescaled pressure and pore space velocity, respectively, and (10)-(11) can be thought of as rescaled Stokes equations to be solved within the -cell. In (10), is the Kronecker delta function ( for , and for ).
Upon further mathematical steps (collecting of like order terms, taking the spatial average over the cell, and applying the divergence theorem and the periodicity of the cell), the spatial average of in the cell is found to be divergence free; . This implies that the macroscopic flow equation is satisfied by Darcy’s equation (1), that is,where is the macroscale pressure and is the macroscale velocity, when the full permeability tensor for an -dimensional problem is defined aswhere is the volume of the -cell (solid and fluid space combined).
In this work, we use the staggered grid finite difference method (i.e., the marker and cell or MAC scheme [36]) for the numerical solution of the rescaled Stokes problem defined by (10)–(12) in the -cell. Our numerical implementation using a matrix-oriented approach is explained in [37]. To illustrate the numerical scheme, let us consider the following Stokes problem (in its original form before being rescaled by homogenization) in two-dimensional space: find the pore space velocity and the pressure such thatwhere the velocity has a unique solution while the pressure is unique up to an additive constant. Here and denote the length of the rectangular domain in the - and -directions, respectively.
To construct a staggered grid finite difference scheme for the Stokes problem, we introduce a possibly nonuniform grid of as follows: The grid is comprised of fluid and solid cells (each cell being fully occupied by one phase), and the computational scheme acts uniquely on the computational cells belonging to the fluid. The pressure at the (fluid) cell center is approximated by . Similarly, we approximate the -velocity at each -edge (fluid) center by , and the -velocity at each -edge (fluid) center by .
By integrating the -component of (15) over cell , we have
We then approximate those edge integrals by midpoint values and the volume integral by its center value and divide both sides by .
By introducing the following forward and backward difference operators:and similar notations for operations in the -direction, we can express the results as
Similarly, by integrating the -component of (15) over cell and then applying midpoint quadrature rule, we can obtain
By integrating (16) over cell , we obtain its discrete form
Eqns. (21)–(23) form a linear equation system for the pressure and velocity, and the equation system defines the MAC scheme for the Stokes problem of (15)–(17). Note that the no-flow boundary condition supplies the zero velocity component that acts normal to any fluid-solid interface, and periodic conditions are applied along the domain’s external boundaries.
We apply this numerical scheme to solve the set of rescaled Stokes equations in Eqns. (10)–(12). Our computer-generated media samples are 2-dimensional (2D); thus the set is solved with different forcings where and . Upon obtaining the rescaled pore space velocity in the -cell, we then compute the full permeability tensor by (14). If our computer-generated media are characterized by a nondiagonal tensor, we determine the principal components of the permeability tensor by diagonalization [11, 38]. The notation used in our results section is and for the permeability tensors in coordinate systems and , respectively, that is,
3.2. Hydraulic Tortuosity
The approach commonly used in literature to compute hydraulic tortuosity is presented in Koponen et al. [19, 26], and further validation of this approach is given in Duda et al. [13]. In this approach, the hydraulic tortuosity or a particular direction is computed based on the pore space fluid velocities driven by a force. For example, the hydraulic tortuosity of flow that is driven in the direction iswhere is the spatial average in the unit cell domain , is the magnitude of the fluid velocities (i.e., the fluid speed), and is the velocity component in the direction. The spatial average of in iswhere is the volume of (pore and solid space combined), is the fluid space, and the integral is taken over the fluid space only. The spatial average for is taken in the same way. The tortuosity in the direction is computed in a similar manner as follows:where the velocity is the solution for flow driven in the direction. We note that (25) and (27) correspond to the hydraulic tortuosity as we have defined it in Section 2.2; that is . If the term which appears in the Kozeny-Carmen equation is labelled as the hydraulic tortuosity, then (25) and (27) should be squared.
While this method has been used in many works (i.e., [26, 33, 39, 40]), other approaches to compute hydraulic tortuosity have been proposed. For example, Matyka et al. [24] took measurements of streamlines that passed through a given cross-sectional area in the domain with a constant flux (i.e., not all streamlines in the domain would be measured). The aim of this approach is to capture the hydraulic tortuosity of the main conducting channels in a pore geometry. Also, Ahmadi et al. [32] formulated analytical functions for hydraulic tortuosity based on volume averaging of mass balance equations.
In our work we follow Koponen et al.’s [19, 26] approach as given by (25)–(27); however we do not compute the pore space velocity fields and explicitly. Instead, during the permeability computation for our 2D pore geometries, we obtain the rescaled pore space velocity fields in the -cell, where and . It can be seen that and are related to and , respectively, for flow that is driven by and that and are related to and , respectively, for flow that is driven by . Thus, it is mathematically equivalent to compute the hydraulic tortuosities bywhere refers to the domain of the -cell, and the spatial averages of are computed in the same way as shown by (26). Since the numerical scheme computes on the fluid edge-centers, we compute the fluid cell-center quantities by averaging across the grid cells, where on the fluid-solid interfaces as per the boundary condition in (12).
As mentioned, we diagonalize the permeability tensor to obtain the principal permeability components and the principal directions and . Thus, in order to make an appropriate comparison between direction-specific permeability and hydraulic tortuosity, we compute the hydraulic tortuosity in these principal directions. For example, in the principal direction iswhere and where , , and are obtained by vector transformations from . in the other principal direction is similarly obtained. The notation we use in our results section is and .
3.3. Effective Diffusion Coefficient and Diffusive Tortuosity
The problem is similarly defined by the two scales illustrated in Figure 1. Within the pore space, the diffusive and convective transport of a species of concentration is given bywhere is divergence free (5), implies the variables are functions of the microscale, and is the molecular diffusion coefficient. The last two equations give the boundary conditions on the fluid-solid interface. Using similar homogenization steps as presented in Section 3.1, the multiscale expansionand (8) are substituted into the governing transport equations (30). Upon collecting like order terms, the second-order term of is where is the macroscale concentration and satisfieswhere (34) gives the boundary condition and is a unit vector normal to the fluid-solid interface (positive in the direction away from the fluid space). Eqns. (33) and (34) correspond to the assumption of an isotropic molecular diffusion coefficient, that is, . Upon further algebra and averaging over the -cell, the macroscale transport equation is obtained as follows:where the effective diffusion coefficient is given byIntuitively, the more tortuous the porous media, the slower the rate of effective diffusion; thus a commonly used relationship between effective and molecular diffusivity iswhere is the diffusive tortuosity [41].
The above approach to compute the effective diffusion coefficient and the diffusive tortuosity has been used in several previous studies, which formulate (33)–(36) using the method of volume averaging [1, 42–45] or by the method of homogenization [27]. Other approaches to compute diffusive tortuosity include simulation of the diffusion process using a random walk [25] or numerically solving the microscale diffusion equation using lattice Boltzmann modeling [46]. Typically, the effective diffusivity is computed and then used to calculate diffusive tortuosity by the inverse relationship shown in (37).
By (36) and (37), diffusive tortuosity is computed bywhere and where is the solution to (33) with the boundary condition in (34). This is the same formulation that was used in Sun et al. [27] by the method of homogenization and in Kim et al. [43] and Valdés-Parada et al. [1] by the method of volume averaging.
The full diffusive tortuosity can also be diagonalized using the same procedure outlined in Section 3.1 to obtain the components and corresponding to the principal directions. The rotation required to obtain the principal components of is not necessarily identical to the rotation required to obtain the principal components of , as will be shown by one of the following pore geometry examples (randomly distributed squares).
4. Results
While real geological formations are typically of low porosity (), we have generated several idealized geometries within the porosity range of . These idealized geometries and porosities have been used in literature (see Table 1) and thus provide ground for comparison and opportunity to make useful comparisons of the effective properties in high porosity representations of porous media.
Permeability, hydraulic tortuosity, and diffusive tortuosity are computed for each pore geometry using the methods outlined in Sections 3.1, 3.2, and 3.3, respectively. Each pore geometry is generated by defining certain input parameters (i.e., radius of solid circle, side length of square, and density of randomly distributed squares), and these inputs are summarized in Table 4. By using a range of values for the input parameters, we obtain unit cells which we treat as representative elementary volumes over a range of porosities. After computing the permeability and tortuosity, we plot the data to observe permeability-porosity and tortuosity-porosity trends.
4.1. In-Line Array of Uniform Shapes
The unit cell, or -cell, in Figure 1(b) illustrates our first type of generated pore geometry, and the shape in the center of the cell is either a circle or square. A range of porosities of the unit cell are obtained by changing the radius of the solid circle or the length of the solid square centered in the cell.
4.1.1. Permeability and Hydraulic Tortuosity
The full permeability tensor was computed for , and results are plotted in Figure 2. Following [14], we computed and plotted the analytical solution from Gebart [47] which is explained in Appendix A. Permeability is isotropic (i.e., ) for this pore geometry. The range of our simulated data fits well against the analytical solution. We computed our permeability data on increasingly finer meshes until the permeability converged within a relative error. To reach this convergence criteria after starting with a mesh size of , the lowest porosity cell required a refined mesh of , and the highest porosity cell required . The shape of the circular edge was refined at each mesh refinement in order to represent the circle as realistically as possible and to minimize artifacts caused by discretization.
We obtained the hydraulic tortuosity in the and directions, for porosities within the range , with the same convergence criteria used to obtain permeability. Figure 3 illustrates the fluid speed and direction for a unit cell of , and the hydraulic tortuosity in both and directions were computed to be . The fluid speed is highest in the pore space that is uninterrupted by the solid, and the flow lines indicate the direction is parallel to the macroscopic flow. The flow lines diverge in the region that is interrupted by the solid; however the fluid speed is essentially zero in that pore space. So, while these flow lines appear to illustrate a tortuous pathway through this geometry, the tortuosity in either direction is in fact calculated to be ≈1 because the highest speeds are acting along a nontortuous pathway. In the pore space that is not restricted or interrupted in the direction parallel to the driving force, the fluid velocity is able to develop into the parabolic profile that characterizes channel flow.
4.1.2. Diffusive Tortuosity
We obtained the full diffusive tortuosity tensor on increasingly refined meshes until the tortuosity converged within a relative error. Starting with a mesh of , the lowest porosity required a refined mesh of , while the highest porosity required to meet the convergence criteria. Solid edges were also refined with mesh refinement, as mentioned previously. In Figure 4, we plot the convergence behavior of for a select number of porosity cases in order to illustrate that each porosity case convergences with similar behavior. Figure 5(a) illustrates the fields of ( or ) for a unit cell of , which were obtained by solving the homogenized problem described in Section 3.3. The resulting diffusive tortuosity tensor iswhere . Notice the tensor is symmetric and diagonal (i.e., to machine precision) and isotropic (i.e., ). To illustrate the passive diffusive transport through this geometry, we have illustrated the diffusion speed and direction in Figure 5(c), which was computed from the concentration field illustrated in Figure 5(b). This concentration field satisfies with Dirichlet boundary conditions on the right and left boundaries such that the concentration gradient across the unit cell is one.
(a)
(b)
(c) Diffusion speed and direction
4.1.3. Comparison of Tortuosities for In-Line Array Unit Cells
The comparison between our isotropic results for and illustrated in Figure 6 allows us to emphasize that these tortuosity types are quantitatively different. This may be surprising at first, since the streamlines in Figures 3 and 5(c) appear to indicate that the flow fields are relatively analogous. However, the fluid speeds shown in Figure 3 are such that the fastest moving flow exists along the nontortuous pathway, while the fastest passive diffusive transport (or highest concentration gradient) shown in Figure 5(c) exists along the circle’s rounded fluid-solid boundary. This difference in speeds (or streamline weights) leads to a quantitative difference between these two tortuosity types.
(a)
(b)
To further make our point, we quantitatively compare the hydraulic and diffusive streamlines in Figures 3 and 5(c). The length of 12 streamlines with starting positions given by is tabulated in Table 5. These lengths are computed by , where is the length of the streamline and for the length of the pore geometry. This table reveals that the length (or tortuosity) of a hydraulic streamline is relatively similar to the diffusive streamline with the same starting location. However, the hydraulic speeds are such that the fastest moving flow exists along the nontortuous pathway, while the fastest diffusive flow (or highest concentration gradient) exists along the circle’s rounded fluid-solid boundary. These speeds mean that individual streamlines contribute to the effective tortuosity with different weights; since the weight of an individual hydraulic streamline is different from the diffusive streamline, the overall hydraulic and diffusive tortuosity of the unit cell are expected to be quantitatively different.
In Figure 6, we plot the approximate analytical solution from Rayleigh [48], which is (see Table 3). Given a system comprised of a periodic array of cylinders, one can imagine that there is a limit as to how small the porosity can become before the pore space is no longer conducting. We find that the diffusive tortuosity data for in-line array of squares fits to this solution with an expected convex shape (comparison between Rayleigh’s [48] analytical solution and simulated data for in-line array of squares was shown in Figure 4 of Valdés-Parada et al. [1]), but the data for in-line array of circles deviates as its critical porosity () is approached. At first glance, one might think the reason for this discrepancy is due to an inadequate mesh refinement: as the solid circle becomes larger, the fluid space between adjacent solids becomes smaller, and a finer mesh may be required to obtain an accurate solution. However, the convergence of our in-line array of circles results was demonstrated in Figure 4. Furthermore, our results match quite closely those obtained by F. J. Valdés-Parada (provided via personal communication, December 3, 2016) using the finite element method. As such, our diffusive tortuosity rather is exhibiting validated convex and concave trends.
In terms of the hydraulic tortuosity, we observe that it is independent of porosity (or a very weak function of ) for this isotropic (symmetrical) geometry. Regardless of the solid circle’s radius, the fluid speeds are highest in the pore space that lies parallel to the driving force direction, as noted in Figure 3. Since the fluid speeds are highest along a nontortuous pathway, the tortuosity of the macroscopic flow through the unit cell is computed to be close to one. We can expect the scenario to be different for a geometry that contains solids which interrupt the main flow path, as seen in the next example of staggered solids.
4.2. Staggered-Array of Uniform Shapes
Figure 7 illustrates our second type of pore geometry. Unlike the previous pore geometry, staggered-array of uniform shapes is an anisotropic geometry. Input parameters to define this geometry are illustrated in Appendix B (Figure 8), and for the following staggered-array examples, and . By varying the solid square length, a range of porosities are obtained.
(a) = 0.36
(b) = 0.64
(c) = 0.84
4.2.1. Permeability and Hydraulic Tortuosity
The permeability for our second pore geometry was computed for a range of porosities. Starting with an initial mesh resolution of , the meshes required refinement up to a maximum of in order to meet the convergence criteria of .
This pore geometry configuration is anisotropic (i.e., ), and the trend of the principal permeabilities over the porosity range tested is shown in Figure 9. The analytical solution for an in-line array of cylinders (i.e., Gebart [47]) is plotted against our staggered-array permeability for comparison, while a fit is not expected. The permeability is characterized by an anisotropic ratio as follows:where for all porosities studied, as shown in Figure 10.
The hydraulic tortuosity in the and directions was computed using the same convergence criteria used to obtain permeability. Our results for a staggered-array unit cell of are shown in Figure 11 and the resulting hydraulic tortuosities were and . In this example, the pore geometry is such that there is not a straight (or uninterrupted) flow path when the flow is driven in the direction. However, there is a straight flow path when the flow is driven in the direction. Due to this characteristic, we observe that the hydraulic flow is more tortuous in the direction than in the direction, which is confirmed by the anisotropic ratio for hydraulic tortuosity plotted in Figure 10. By observing Figure 11, we notice that the fluid speeds are highest in the pore regions located between two square corners. We also notice in Figure 11(a) that the fluid speed is close to being constant along a sinusoidal flow path which appears to occupy most of the pore space, except towards the fluid-solid boundaries. The horizontally driven flow is interrupted by the center square and causes the flow path to diverge. However, in Figure 11(b), the vertically driven flow is not interrupted to the extent that would cause the flow path to become sinusoidal. Instead, the geometry of the pore space allows for a straight flow path to develop. We do not notice this sinusoidal flow path of constant speed, rather we notice that the speed is more localized (i.e., speed is highest between adjacent solid corners and lowest along solid square sides).
(a) Flow in direction
(b) Flow in direction
4.2.2. Diffusive Tortuosity
We obtained the full diffusive tortuosity tensor for this staggered geometry using a convergence criteria of . Starting with a uniform mesh of , the lowest porosity required a refined mesh of while the highest porosity required in order to meet the convergence criteria. Results for a staggered-array unit cell of are shown in Figure 12, and the diffusive tortuosity tensor iswhere the tensor is symmetric and diagonal (i.e., to machine precision) and anisotropic (i.e., ). In this geometry of staggered squares, , and thus passive diffusive transport is more tortuous in the direction due to an interrupted flow path.
(a)
(b)
(c) Diffusion speed and direction
(d)
(e)
(f) Diffusion speed and direction
4.2.3. Comparison of Tortuosities for Staggered-Array Unit Cells
Our results for hydraulic and diffusive tortuosity over the range of porosities are plotted in Figure 13. Diffusive tortuosity data from Ryan et al. [49] and Kim et al. [43] (both obtained from Kim et al. [43]) are also plotted for comparison. While the works of Ryan et al. [49] and Kim et al. [43] were based on using the method of volume averaging to obtain the boundary-value problem for diffusivity presented in Section 3.3, we are able to compare our results to theirs since both homogenization and volume averaging approaches lead to the same boundary-value problem. To solve the boundary-value problem presented in Section 3.3, Ryan et al. [49] employed the finite difference method while Kim et al. [43] employed the boundary element method. We note that one of Kim et al.’s [43] data points shown in Figure 13(b) unexpectedly corresponds to a diffusive tortuosity of less than 1. After confirming this value from Table 6 of their work (where , /), we do not have an explanation as to how could be smaller than 1, other than a possible typographical error. Since they reported values for the parameter rather than explicitly, it may not have been easily apparent to them that one of their values corresponded to . (Refer to Appendix B where we explain more about the notation as well as compare our numerical results against more of their results.) Regarding the diffusive tortuosity, our simulated data fits closely to these two literature data; however, in Figure 13(a) we observe our data diverges at lower porosities (). While the exact reason for this discrepancy is not clear, possible causes may be attributed to the difference in numerical methods employed or the geometry of the staggered-array (especially at lower porosities that contain very narrow flow channels). Figure 13(b) shows a closer fit between our data to the literature data, in comparison with the fit shown in Figure 13(a). Since the narrowest channels exist at the lowest porosity and are oriented in the direction, this observation implies our diffusive tortuosity data may be sensitive to the mesh resolution between adjacent solid boundaries. Despite the difference between our data to that reported in literature, the anisotropic nature of this geometry is evident from our results. Over the range of porosities, we observe that the hydraulic flow is predominantly parallel to the driving force, unless an obstacle prevents a straight flow path. Thus, for all porosities as shown in Figure 13, and is a very weak function of porosity. The anisotropic ratios of the tortuosities are and are plotted in Figure 10, along with the permeability ratio previously defined. Over the whole porosity range, both anisotropic ratios of the tortuosities are less than one, indicating the hydraulic flow and passive diffusive transport are more tortuous in the direction than in the direction. The hydraulic tortuosity shows a greater degree of anisotropy than diffusive tortuosity because for all porosities while increases towards the lower porosities.
(a) Tortuosity in horizontal direction
(b) Tortuosity in vertical direction
4.3. Randomly Distributed Squares
Our third pore geometry generated comprised of randomly distributed squares, which is a type commonly used in literature (e.g., [13, 19, 26]). In those works, squares are freely overlapping; however, in our work, our square size is equivalent to the cell size of the initial mesh resolution, and thus no squares overlap. In this geometry, periodic boundary conditions were implemented by checking for the existence of any solid-fluid interfaces along the external boundaries of the domain, and, if detected, no-flow boundary conditions were applied at the interface (i.e., the velocity component acting normal to the fluid-solid interface was assigned a value of zero in the system of equations). We imposed the constraint given by Koza et al. [34] (i.e., the ratio of obstacle length to domain length < 0.01) which was recommended to avoid anisotropic-related statistical errors. This same constraint was also respected in Duda et al. [13]. Our geometries ranged from (effective) porosities of . Example pore geometries are shown in Figure 14. Our algorithm to generate these pore geometries included a step to check for any fluid sites that may be isolated from the flow path, which if found were changed to be solid sites. By filling in the isolated fluid sites, we were able to represent the effective porosity of the geometry, as noted in Koponen et al. [26]. Since lower porosity geometries contained more solid sites than higher porosities, they were especially susceptible to isolated fluid sites. Filling in the isolated sites created large solids within the pore geometry (especially noted in Figure 14(a)), and the impact of these solids on results will be discussed. Due to the nature of this randomly generated geometry, we used between 8 and 11 realizations of a given porosity between 0.46 and 0.95 and 5 realizations for the porosity of 0.97. Thus the following results are presented as averaged quantities with error bars to indicate the spread between the maximum and minimum quantities computed.
(a)
(b)
(c)
(d)
(e)
4.3.1. Permeability and Hydraulic Tortuosity
We computed permeability using a convergence criteria of . Starting with an initial mesh of , the lowest porosity required a refined mesh of and the highest porosity required to meet the convergence criteria.
The full permeability tensors (in coordinate system ) for these pore geometries were symmetric and nondiagonal; thus a transformation was done to obtain the diagonal tensor in the principal coordinate system . The full permeability components and diagonal permeability components are plotted in Figure 15. The nature of the anisotropic ratio of the diagonal permeability tensor will be discussed and shown later in Figure 16(b).
(a) Full tensor components, in
(b) Diagonalized tensor components, in
(a) Rotation angle,
(b) Anisotropic ratios
We computed hydraulic tortuosity from the flow field results which were obtained using the convergence criteria for permeability previously mentioned. For a configuration case of , results for hydraulic flow are shown in Figure 17 and the computed hydraulic tortuosities are and . From Figure 17, the tortuous nature of the fluid pathways is evident. We note that a few “hot spots” of fastest moving fluid (with speeds ) exist in the geometry, while the rest of the pore space is characterized by moderate () to slow () moving fluid. The high speeds exist within the pore geometry’s narrow channels, and, upon magnifying these regions, we observed a parabolic velocity profile indicative of channel flow. While Figure 17 illustrates fluid flow in the direction only, similar characteristics were found for flow in the direction; thus we have omitted those illustrations. The computed hydraulic tortuosity at this porosity is essentially isotropic; that is, . So while the geometry is comprised of randomly distributed squares, the overall geometry exhibits an isotropic hydraulic tortuosity. (However, this is not the case at lower porosity configurations, as the anisotropic ratios which we will discuss later using Figure 16(b) will indicate.) In agreement with results presented in [13], the trend shown in Figure 18 indicates that the flow pathways are more tortuous as the geometry is comprised of more solid squares.
(a) Fluid speed
(b) Close-up of boxed region in (a)
4.3.2. Fit to Kozeny-Carman Equation
As presented in Section 2.2, the Kozeny-Carman equation relates permeability to porosity , hydraulic tortuosity , specific surface area , and a fitting coefficient (which we refer to as the shape factor). Specific surface area for the range of porosities was computed by (3) and is shown in Figure 19. The specific surface area increases as the porosity decreases because the geometry is comprised of more solid squares. However, once the porosity is lowered to the point where isolated fluid sites must be converted to large solid areas, the fluid-solid interfacial area is reduced along with the specific surface area. We computed the shape factor by (2) for this particular pore geometry, using our permeability, hydraulic tortuosity, and specific surface area data presented in Figures 15(a), 18, and 19, respectively. The trend for over the porosities is plotted in Figure 20, and we observe that the shape factor does not change significantly within a porosity interval of , which was similarly reported in Duda et al. [13].
4.3.3. Diffusive Tortuosity
We obtained diffusive tortuosity for this geometry using different convergence criteria for three ranges of porosities due to high mesh refinement requirements. Convergence criteria of , , and were used for porosity ranges of , , and , respectively. The geometries were generated using a mesh of , and the convergence criteria required a refinement of for the lowest porosity and for the highest porosity. Results for a geometry of are shown in Figure 21, and the resulting diffusive tortuosity tensor isNotice the tensor is symmetric (i.e., ) but not diagonal (i.e., the off-diagonal components are nonzero). As such, we diagonalized the full diffusive tortuosity tensor and results for the whole porosity range studied are presented in Figure 22.
(a)
(b)
(c)
(d) Close-up of boxed region in (c)
(a) Full tensor components, in
(b) Diagonalized tensor components, in
4.3.4. Comparison of Tortuosities for Randomly Distributed Squares
Our hydraulic and diffusive tortuosity data are plotted in Figure 23 for a range of porosity values, along with some of the trends from literature that we reported in Tables 2 and 3. The trend from Koponen et al. [26] was obtained for hydraulic tortuosity for randomly distributed squares. The trend from Mackie and Meares [50] was for diffusive transport of electrolytes through a membrane, and their expression for the tortuosity-porosity trend was credited in the past as one of the most successfully employed correlations for membrane systems [51]. The trend from Weissberg [52] was also for diffusive transport, through a geometry comprised of uniform spheres. It was given as an upper bound for effective diffusion and thus could be considered as the lower bound for diffusive tortuosity since they are inversely related. Our hydraulic data fits moderately close with the Koponen et al.’s [26] trend. Our diffusion data follows a similar trend to Mackie and Meares [50] until and remains above Weissberg’s [52] lower bound. In general, our data demonstrates that diffusive tortuosity is not equivalent to hydraulic tortuosity and can be up to ten times greater, especially at low porosities.
(a) Tortuosity-porosity trends
(b) Close-up of data
4.3.5. Induced Anisotropic Properties at Lower Porosity Geometries
Due to the method we used to generate the randomly distributed squares geometry, anisotropic properties were introduced into the pore structure, especially at lower porosities. This is evident in Figure 24, where low and medium porosities are compared. The lower porosity has large, nonuniform solid shapes within the geometry. To represent the connected pore space or effective pore space [26], we filled in the isolated fluid sites, which resulted from the random distribution of squares in the domain. The result is a nonuniform distribution of solids in the geometry and an induced degree of anisotropy that is shown in Figure 16(b). (Note that Figure 16(a) shows that the rotation required to obtain the diagonal tensors, and , is independent of porosity.) At the higher porosities, the degree of anisotropy of all three properties is moderately close to one, and thus we could consider those geometries as exhibiting isotropic behavior. As the porosity is lowered, the permeability and diffusive tortuosity become anisotropic and the degree of anisotropy of the permeability is higher than that of the diffusive tortuosity. The degree of anisotropy of hydraulic tortuosity is less than diffusive tortuosity and is close to being isotropic for most of the porosity range. A possible reason for this finding is that hydraulic flow develops preferential pathways where the fluid speed is highest along a pathway that is unhindered by solids, as shown in Figure 25(b). If vertical and horizontal preferential pathways are geometrically similar, an isotropic hydraulic tortuosity is exhibited. On the other hand, the diffusive speeds are highest adjacent to the fluid-solid boundaries, as seen in Figure 25(d); thus the geometry of the pore structure appears to impact the diffusive tortuosity more than it impacts the hydraulic tortuosity. (This point is also illustrated by comparing the flows in a higher porosity, in Figures 17(b) and 21(d).)
(a)
(b)
(a) Fluid speed
(b) Close-up of boxed region in (a)
(c) Diffusion speed
(d) Close-up of boxed region in (c)
5. Conclusion
In this study, we computed the effective properties known as permeability, hydraulic tortuosity, and diffusive tortuosity, using well-known and commonly employed formulations based on the method of homogenization. While we noted that past work has recognized the difference between tortuosity types, we emphasized that a few studies have failed to distinguish the difference and have mistakenly used them interchangeably. In this work, hydraulic tortuosity was computed based on the approach given in Koponen et al. [19, 26]. Their approach uses pore space velocity fields ; however we used the rescaled velocity fields which we obtained by solving a set of rescaled Stokes equations formulated by homogenization ((10)–(12)). Using in (28) is mathematically equivalent to using in (25) and (27). Diffusive tortuosity was computed by solving the homogenized or boundary-value problem given for the effective diffusion coefficient, as employed in Kim et al. [43], Valdés-Parada et al. [1], and Sun et al. [27].
We generated several different pore geometries, including both in-line array and staggered-array of uniform shapes, and randomly distributed squares. For the in-line array of circles geometry, our simulated permeability data was validated against the analytical solution from Gebart [47], and our simulated diffusive tortuosity data was validated against Rayleigh’s [48] trend. We studied the anisotropy of the computed properties and fit data from one of the pore geometries to the Kozeny-Carman equation to obtain the shape factor .
The main findings from the geometries studied are as follows:(1)Hydraulic tortuosity is not equal to diffusive tortuosity in the same pore geometry. In the in-line array of uniform shapes (either circle or square), hydraulic tortuosity was weakly dependent on porosity, while diffusive tortuosity was almost linearly related to porosity (recall Figure 6). In the randomly distributed squares geometry, the diffusive tortuosity was greater (up to a factor of ten) than hydraulic tortuosity as porosity decreased (recall Figure 23). In all examples, hydraulic speeds were highest along the mid-channel space between adjacent solids and formed a parabolic velocity profile, while diffusive speeds were highest along the fluid-solid boundary which could be extremely irregular or complex (as seen in the randomly distributed squares geometry); we suspected this boundary impact on flow to be the main reason why these two tortuosity types were quantitatively different in the same pore geometry.(2)Results from the three different pore geometry configurations used in this work indicate that for the majority of the porosity range tested. However, we do not claim that this relationship is necessarily true for all types of porous media. For example, a pore network model that was considered in Ghanbarian et al. [8] led to the conclusion that . Thus we simply note an equality held true for the three types of synthetic porous media tested in our work.(3)Related to the first point, hydraulic and diffusive tortuosity (as well as permeability) can exhibit different anisotropic behavior in the same pore geometry. In the staggered-array of uniform shapes (square), we observed that the hydraulic tortuosity had a greater degree of anisotropy compared to diffusive tortuosity’s anisotropy (recall Figure 10). This behavior was due to the nature of the geometry, which allowed for a nontortuous hydraulic flow in one of the principal directions. In the randomly distributed squares, the large nonuniform shapes that were introduced into the geometry as a result of isolated fluid sites led to an induced anisotropic behavior; at higher porosities, all properties exhibited essentially isotropic behavior, but at lower porosities, and became anisotropic by varying degrees. It is interesting to note that displays the least degree of anisotropic behavior; this finding may be related to the existence of preferential pathways or regions of high-speed flow through minimally tortuous pathways within the pore geometry.(4)Two geometries (i.e., staggered-array of squares, randomly distributed squares) demonstrated that when the permeability is anisotropic, the effective diffusion coefficient (which is inversely proportional to diffusive tortuosity) can also be anisotropic however not necessarily by the same degree. In general, the degree of anisotropic permeability was greater than the degree of anisotropic diffusion in the same pore geometry (recall Figures 10 and 16(b)). This finding has important implications for flow and transport modeling at reservoir-scale because the use of Darcy’s equation and a transport equation (with a diffusive flux term) requires permeability and the effective diffusion coefficient as input parameters, respectively.(5)A few qualitative statements can be made regarding the flow behavior demonstrated in some of the examples (i.e., staggered-array of squares and randomly distributed squares). Hydraulic tortuosity and permeability are generally related; the more tortuous the flow pathway is in a direction, the slower the fluid speed is and the less permeable it is in that direction. Lower porosity geometries are characterized by slower fluid speeds and more tortuous hydraulic pathways, in addition to more tortuous diffusive pathways. By (37), which shows the effective diffusion coefficient is inversely proportional to diffusive tortuosity, we can deduce that the more tortuous the diffusing flow in a direction, the slower the rate of (effective) diffusive transport.
Appendix
A. Validation of Permeability Calculation
To validate our numerical implementation and computation of permeability, we compared our results with the analytical solution given by Gebart [47]. This analytical solution of permeability corresponds to an in-line array of infinitely long cylinders. The porosity is computed as a function of the cylindrical radius as follows:where are the fluid and solid volumes, respectively, is the total volume, and the unit cell is . Thus by varying , the porosity is also varied. Gebart’s [47] analytical solution for permeability (normalized by ) as a function of porosity iswhere is the geometric factor which depends on the fiber arrangement, and is the critical porosity (percolation threshold). For a quadratic fiber arrangement, the geometric factor is , and the percolation threshold is for an array of cylinders.
B. Validation of Diffusive Tortuosity Calculation
For benchmarking purposes, we computed the diffusive tortuosity for staggered unit cells and compared our results to those reported by Kim et al. [43], who used the boundary element method (BEM). In our work, we used finite difference. The method to compute the diffusive tortuosity tensor was presented Section 3.3, and for a 2D problem the tensor components are computed bywhere . Note that diffusive tortuosity is the inverse of the diffusive tortuosity factor, that is, . Recalling (37), the diffusive tortuosity factor is related to the effective diffusion coefficient by where . The form of (B.2) helps to make appropriate comparisons in Table 6. The input parameters required to construct the staggered unit cell are porosity , ratio, and ratio, and it is constructed to meet the configuration presented in Figure 8. Kim et al. [43] performed measurements on many staggered unit cell cases, and Table 6 presents a comprehensive comparison of our results against their data.
Disclosure
Rebecca Allen is now at SINTEF Digital, Mathematics and Cybernetics, Forskningsveien 1, 0314 Oslo, Norway.
Competing Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
The authors gratefully acknowledge that the research reported in this publication was supported by funding from King Abdullah University of Science and Technology (KAUST). The authors gratefully acknowledge and thank Dr. Francisco J. Valdés-Parada for providing the diffusive tortuosity results for in-line array of circle and square geometries, which were used for comparison in Figure 6 with permission.