Research Article  Open Access
Computing and Comparing Effective Properties for Flow and Transport in ComputerGenerated Porous Media
Abstract
We compute effective properties (i.e., permeability, hydraulic tortuosity, and diffusive tortuosity) of three different digital porous media samples, including inline array of uniform shapes, staggeredarray 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 reservoirscale 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 porescale 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 wellknown and important properties used in Darcy and reservoirscale 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].
Porescale 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 computergenerated porous media samples; real refers to a digital sample obtained by imaging and conversion to binary form, and computergenerated 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 porescale simulation to obtain effective properties dates back to Cancelliere et al. [7] and possibly earlier.

In regard to obtaining effective properties at the porescale, 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 KozenyCarman 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 (inline array of uniform shapes, staggeredarray 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 porescale (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 porescale 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 opensource 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 KozenyCarman (KC) equation [17, 18].
2.2. KozenyCarman Equation
An idealized pore geometry of parallel cylindrical channels was used to formulate the KozenyCarman (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 crosssectional 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 fluidsolid interface area to the total volume in the domain:Notice that (2) is not written as directionspecific. In fact, Carrier III [20] mentions that a limitation of the KC equation is that it does not explicitly account for anisotropy (a directionspecific 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 KozenyCarman equation could still apply to directionspecific 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 platelike obstacles and used the diffusive tortuosity results to compute fitting coefficients present in Koponen et al.’s [26] tortuosityporosity 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ésParada et al. [1] stated that its quantification becomes more complex when more than just passive diffusion is occurring in the system. ValdésParada 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 TortuosityPorosity 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.

 
referenced in [57]. referenced in [30, 31]. 
In the same way that the KC equation is specific to a given class of porous media, any porositytortuosity 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 platelike 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 Xray 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ésParada 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 tortuosityporosity 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 directionalspecific 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 fluidsolid interface is given bywhich comes from the noflow and noslip boundary conditions. Although the slip boundary condition is necessary in certain contexts (such as gas flow and nonNewtonian fluid), the noslip boundary condition is a valid assumption for viscous flow in porous media and is widely used in modeling studies of porescale 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 secondorder pressure term and the thirdorder velocity term aresuch that and satisfywith the fluidsolid 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 matrixoriented 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 twodimensional 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 noflow boundary condition supplies the zero velocity component that acts normal to any fluidsolid 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 computergenerated media samples are 2dimensional (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 computergenerated 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 KozenyCarmen 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 crosssectional 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 edgecenters, we compute the fluid cellcenter quantities by averaging across the grid cells, where on the fluidsolid 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 directionspecific 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 fluidsolid 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 secondorder term of is where is the macroscale concentration and satisfieswhere (34) gives the boundary condition and is a unit vector normal to the fluidsolid 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ésParada 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 permeabilityporosity and tortuosityporosity trends.

4.1. InLine 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 InLine 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 fluidsolid 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 fluidsolid 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 inline array of squares fits to this solution with an expected convex shape (comparison between Rayleigh’s [48] analytical solution and simulated data for inline array of squares was shown in Figure 4 of ValdésParada et al. [1]), but the data for inline 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 inline array of circles results was demonstrated in Figure 4. Furthermore, our results match quite closely those obtained by F. J. ValdésParada (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. StaggeredArray of Uniform Shapes
Figure 7 illustrates our second type of pore geometry. Unlike the previous pore geometry, staggeredarray of uniform shapes is an anisotropic geometry. Input parameters to define this geometry are illustrated in Appendix B (Figure 8), and for the following staggeredarray 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 inline array of cylinders (i.e., Gebart [47]) is plotted against our staggeredarray 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 staggeredarray 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 fluidsolid 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 staggeredarray 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 StaggeredArray 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 boundaryvalue 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 boundaryvalue problem. To solve the boundaryvalue 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 staggeredarray (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 solidfluid interfaces along the external boundaries of the domain, and, if detected, noflow boundary conditions were applied at the interface (i.e., the velocity component acting normal to the fluidsolid 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 anisotropicrelated 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) Closeup of boxed region in (a)
4.3.2. Fit to KozenyCarman Equation
As presented in Section 2.2, the KozenyCarman 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 fluidsolid 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 offdiagonal 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) Closeup 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 tortuosityporosity 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) Tortuosityporosity trends
(b) Closeup 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 fluidsolid 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) Closeup of boxed region in (a)
(c) Diffusion speed
(d) Closeup of boxed region in (c)
5. Conclusion
In this study, we computed the effective properties known as permeability, hydraulic tortuosity, and diffusive tortuosity, using wellknown 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 boundaryvalue problem given for the effective diffusion coefficient, as employed in Kim et al. [43], ValdésParada et al. [1], and Sun et al. [27].
We generated several different pore geometries, including both inline array and staggeredarray of uniform shapes, and randomly distributed squares. For the inline 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 KozenyCarman 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 inline 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 midchannel space between adjacent solids and formed a parabolic velocity profile, while diffusive speeds were highest along the fluidsolid 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 staggeredarray 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 highspeed flow through minimally tortuous pathways within the pore geometry.(4)Two geometries (i.e., staggeredarray 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 reservoirscale 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., staggeredarray 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 inline 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ésParada for providing the diffusive tortuosity results for inline array of circle and square geometries, which were used for comparison in Figure 6 with permission.
References
 F. J. ValdésParada, M. L. Porter, and B. D. Wood, “The role of tortuosity in upscaling,” Transport in Porous Media, vol. 88, no. 1, pp. 1–30, 2011. View at: Publisher Site  Google Scholar  MathSciNet
 M. A. Knackstedt, S. Latham, M. Madadi, A. Sheppard, T. Varslot, and C. Arns, “Digital rock physics: 3D imaging of core material and correlations to acoustic and flow properties,” Leading Edge, vol. 28, no. 1, pp. 28–33, 2009. View at: Publisher Site  Google Scholar
 H. Andrä, N. Combaret, J. Dvorkin et al., “Digital rock physics benchmarks—part I: imaging and segmentation,” Computers & Geosciences, vol. 50, pp. 25–32, 2013. View at: Publisher Site  Google Scholar
 S. Torquato and Y. Jiao, “Robust algorithm to generate a diverse class of dense disordered and ordered sphere packings via linear programming,” Physical Review E. Statistical, Nonlinear, and Soft Matter Physics, vol. 82, no. 6, article 061302, 2010. View at: Publisher Site  Google Scholar  MathSciNet
 C. L. Y. Yeong and S. Torquato, “Reconstructing random media,” Physical Review E. Statistical, Nonlinear, and Soft Matter Physics, vol. 57, no. 1, pp. 495–506, 1998. View at: Publisher Site  Google Scholar  MathSciNet
 A. Yang, C. T. Miller, and L. D. Turcoliver, “Simulation of correlated and uncorrelated packing of random size spheres,” Physical Review E  Statistical Physics, Plasmas, Fluids, and Related Interdisciplinary Topics, vol. 53, no. 2, article 1516, 1996. View at: Google Scholar
 A. Cancelliere, C. Chang, E. Foti, D. H. Rothman, and S. Succi, “The permeability of a random medium: comparison of simulation with theory,” Physics of Fluids A: Fluid Dynamics, vol. 2, no. 12, pp. 2085–2088, 1990. 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, pp. 1461–1477, 2013. View at: Publisher Site  Google Scholar
 X. Zhang and M. A. Knackstedt, “Direct simulation of electrical and hydraulic tortuosity in porous solids,” Geophysical Research Letters, vol. 22, no. 17, pp. 2333–2336, 1995. View at: Publisher Site  Google Scholar
 E. SanchezPalencia, “Homogenization method for the study of composite media,” in Asymptotic Analysis II, pp. 192–214, Springer, 1983. View at: Google Scholar
 J. Bear and A. H. Cheng, Modeling Groundwater Flow and Contaminant Transport, vol. 23 of Theory and Applications of Transport in Porous Media, Springer, 2010.
 P. Mostaghimi, M. J. Blunt, and B. Bijeljic, “Computations of absolute permeability on microCT image,” Mathematical Geosciences, vol. 45, no. 1, pp. 103–125, 2013. View at: Publisher Site  Google Scholar  MathSciNet
 A. Duda, Z. Koza, and M. Matyka, “Hydraulic tortuosity in arbitrary porous media flow,” Physical Review E, vol. 84, no. 3, Article ID 036319, pp. 1–8, 2011. View at: Publisher Site  Google Scholar
 A. Ebrahimi Khabbazi, J. S. Ellis, and A. Bazylak, “Developing a new form of the KozenyCarman parameter for structured porous media through latticeBoltzmann modeling,” Computers & Fluids, vol. 75, pp. 35–41, 2013. 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
 A. Nabovati, J. Hinebaugh, A. Bazylak, and C. H. Amon, “Effect of porosity heterogeneity on the permeability and tortuosity of gas diffusion layers in polymer electrolyte membrane fuel cells,” Journal of Power Sources, vol. 248, pp. 83–90, 2014. View at: Publisher Site  Google Scholar
 P. Carman, “Fluid flow through granular beds,” TransactionsInstitution of Chemical Engineers, vol. 15, 1937. View at: Google Scholar
 P. C. Carman, “Permeability of saturated sands, soils and clays,” The Journal of Agricultural Science, vol. 29, no. 2, pp. 262–273, 1939. View at: Publisher Site  Google Scholar
 A. Koponen, M. Kataja, and J. Timonen, “Tortuous flow in porous media,” Physical Review E, vol. 54, no. 1, pp. 406–410, 1996. View at: Google Scholar
 W. D. Carrier III, “Goodbye, Hazen; hello, KozenyCarman,” Journal of Geotechnical and Geoenvironmental Engineering, vol. 129, no. 11, pp. 1054–1056, 2003. View at: Publisher Site  Google Scholar
 J. D. Hyman, P. K. Smolarkiewicz, and C. Larrabee Winter, “Pedotransfer functions for permeability: a computational study at pore scales,” Water Resources Research, vol. 49, no. 4, pp. 2080–2092, 2013. View at: Publisher Site  Google Scholar
 F. D. E. Latief and U. Fauzi, “KozenyCarman and empirical formula for the permeability of computer rock models,” International Journal of Rock Mechanics and Mining Sciences, vol. 50, pp. 117–123, 2012. View at: Publisher Site  Google Scholar
 P. Xu and B. Yu, “Developing a new form of permeability and Kozeny–Carman constant for homogeneous porous media by means of fractal geometry,” Advances in Water Resources, vol. 31, no. 1, pp. 74–81, 2008. View at: Publisher Site  Google Scholar
 M. Matyka, A. Khalili, and Z. Koza, “Tortuosityporosity relation in porous media flow,” Physical Review E, vol. 78, no. 2, Article ID 026306, 2008. View at: Publisher Site  Google Scholar
 T. Ohkubo, “Tortuosity based on anisotropic diffusion process in structured platelike obstacles by Monte Carlo simulation,” Transport in Porous Media, vol. 72, no. 3, pp. 339–350, 2008. 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: Google Scholar
 Z. Sun, X. Tang, and G. Cheng, “Numerical simulation for tortuosity of porous media,” Microporous and Mesoporous Materials, vol. 173, pp. 37–42, 2013. View at: Publisher Site  Google Scholar
 C. David, “Geometry of flow paths for fluid transport in rocks,” Journal of Geophysical Research, vol. 98, no. 7, pp. 12–278, 1993. View at: Google Scholar
 M. Quintard, L. Bletzacker, D. Chenu, and S. Whitaker, “Nonlinear, multicomponent, mass transport in porous media,” Chemical Engineering Science, vol. 61, no. 8, pp. 2643–2669, 2006. 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
 B. P. Boudreau, “The diffusive tortuosity of finegrained unlithified sediments,” Geochimica et Cosmochimica Acta, vol. 60, no. 16, pp. 3139–3142, 1996. View at: Publisher Site  Google Scholar
 M. M. Ahmadi, S. Mohammadi, and A. N. Hayati, “Analytical derivation of tortuosity and permeability of monosized spheres: a volume averaging approach,” Physical Review E  Statistical, Nonlinear, and Soft Matter Physics, vol. 83, no. 2, article 026312, 2011. View at: Publisher Site  Google Scholar
 A. Nabovati and A. C. M. Sousa, “Fluid flow simulation in random porous media at pore level using the lattice Boltzmann method,” Journal of Engineering Science and Technology, vol. 2, no. 3, pp. 226–237, 2007. View at: Google Scholar
 Z. Koza, M. Matyka, and A. Khalili, “Finitesize anisotropy in statistically uniform porous media,” Physical Review E, vol. 79, no. 6, Article ID 066306, pp. 1–7, 2009. View at: Publisher Site  Google Scholar
 Y. Liu and P. K. Kitanidis, “Tortuosity and Archie's law,” in Advances in Hydrogeology, P. K. Mishra and K. L. Kuhlman, Eds., chapter 6, pp. 115–126, Springer, New York, NY, USA, 2013. View at: Google Scholar
 F. H. Harlow and J. E. Welch, “Numerical calculation of timedependent viscous incompressible flow of fluid with free surface,” The Physics of Fluids, vol. 8, no. 12, pp. 2182–2189, 1965. View at: Publisher Site  Google Scholar
 S. Sun, A. Salama, and M. F. El Amin, “Matrixoriented implementation for the numerical solution of the partial differential equations governing flows and transport in porous media,” Computers and Fluids, vol. 68, pp. 38–46, 2012. View at: Publisher Site  Google Scholar  MathSciNet
 J. R. Fanchi, “Mathematics of fluid flow,” in Petroleum Engineering Handbook, L. W. Lake, Ed., chapter 2, pp. 45–76, Society of Petroleum Engineers, 2007, http://petrowiki.org/PEH%3AMathematics of Fluid Flow. View at: 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  Zentralblatt MATH
 M. Icardi, G. Boccardo, D. L. Marchisio, T. Tosco, and R. Sethi, “Porescale simulation of fluid flow and solute dispersion in threedimensional porous media,” Physical Review E—Statistical, Nonlinear, and Soft Matter Physics, vol. 90, no. 1, Article ID 013032, 2014. View at: Publisher Site  Google Scholar
 M. M. Clark, Transport Modeling for Environmental Engineers and Scientists, John Wiley & Sons, 2nd edition, 2009.
 J. A. Ochoa, S. Whitaker, and P. Stroeve, “Determination of cell membrane permeability in concentrated cell ensembles,” Biophysical Journal, vol. 52, no. 5, pp. 763–774, 1987. View at: Publisher Site  Google Scholar
 J. Kim, J. Ochoa, and S. Whitaker, “Diffusion in anisotropic porous media,” Transport in Porous Media, vol. 2, no. 4, pp. 327–356, 1987. View at: Publisher Site  Google Scholar
 A. E. Sáez, J. C. Perfetti, and I. Rusinek, “Prediction of effective diffusivities in porous media using spatially periodic models,” Transport in Porous Media, vol. 6, no. 2, pp. 143–157, 1991. View at: Publisher Site  Google Scholar
 S. Beyhaghi and K. M. Pillai, “Estimation of tortuosity and effective diffusivity tensors using closure formulation in a sintered polymer wick during transport of a nondilute, multicomponent liquid mixture,” Special Topics and Reviews in Porous Media, vol. 2, no. 4, pp. 267–282, 2011. View at: Publisher Site  Google Scholar
 Y. Yong, X. Lou, S. Li, C. Yang, and X. Yin, “Direct simulation of the influence of the pore structure on the diffusion process in porous media,” Computers & Mathematics with Applications, vol. 67, no. 2, pp. 412–423, 2014. View at: Publisher Site  Google Scholar
 B. R. Gebart, “Permeability of unidirectional reinforcements for RTM,” Journal of Composite Materials, vol. 26, no. 8, pp. 1100–1133, 1992. View at: Publisher Site  Google Scholar
 L. Rayleigh, “LVI. On the influence of obstacles arranged in rectangular order upon the properties of a medium,” Philosophical Magazine Series 5, vol. 34, no. 211, Article ID 481502, pp. 481–502, 1892. View at: Google Scholar
 D. Ryan, R. Carbonell, and S. Whitaker, “A theory of diffusion and reaction in porous media,” in Proceedings of the American Institute of Chemical Engineers, vol. 71 of AIChE Symposium Series 202, pp. 46–62, January 1981. View at: Google Scholar
 J. S. Mackie and P. Meares, “The diffusion of electrolytes in a cationexchange resin membrane. I. Theoretical,” Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, vol. 232, no. 1191, pp. 498–509, 1955. View at: Publisher Site  Google Scholar
 S. B. Iversen, V. K. Bhatia, K. DamJohansen, and G. Jonsson, “Characterization of microporous membranes for use in membrane contactors,” Journal of Membrane Science, vol. 130, no. 12, pp. 205–217, 1997. View at: Publisher Site  Google Scholar
 H. L. Weissberg, “Effective diffusion coefficient in porous media,” Journal of Applied Physics, vol. 34, no. 9, pp. 2636–2639, 1963. View at: Publisher Site  Google Scholar
 M. C. Sukop, H. Huang, P. F. Alvarez, E. A. Variano, and K. J. Cunningham, “Evaluation of permeability and nonDarcy flow in vuggy macroporous limestone aquifer samples with lattice Boltzmann methods,” Water Resources Research, vol. 49, no. 1, pp. 216–230, 2013. View at: Publisher Site  Google Scholar
 J. Comiti and M. Renaud, “A new model for determining mean structure parameters of fixed beds from pressure drop measurements: application to beds packed with parallelepipedal particles,” Chemical Engineering Science, vol. 44, no. 7, pp. 1539–1545, 1989. View at: Publisher Site  Google Scholar
 E. du Plessis, S. Woudberg, and J. Prieur du Plessis, “Porescale modelling of diffusion in unconsolidated porous structures,” Chemical Engineering Science, vol. 65, no. 8, pp. 2541–2551, 2010. View at: Publisher Site  Google Scholar
 L. Pisani, “Simple expression for the tortuosity of porous media,” Transport in Porous Media, vol. 88, no. 2, pp. 193–203, 2011. View at: Publisher Site  Google Scholar  MathSciNet
 J. A. OchoaTapia, P. Stroeve, and S. Whitaker, “Diffusive transport in twophase media: spatially periodic models and maxwell's theory for isotropic and anisotropic systems,” Chemical Engineering Science, vol. 49, no. 5, pp. 709–726, 1994. View at: Publisher Site  Google Scholar
Copyright
Copyright © 2017 Rebecca Allen and Shuyu Sun. 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.