Computing and Comparing Effective Properties for Flow and Transport in Computer-Generated Porous Media

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.


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 Darcyand 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][5][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 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.

Study
Real sample Computer-generated sample Porosity range Effective properties computed [16] 3D fibrous material (membrane) 0.45-0.9permeability, hydraulic tortuosity [15] 3D samples: Fontainbleau (F), Berea (B), Carbonate (C) Sphere pack F: 0.147, B: 0.184, C: 0.247, Spheres: 0.343 Permeability, electrical resistivity, among others [14] 2D staggered cylinders, 3D body-centered-cubic spheres 0.25-0.98 Permeability, hydraulic tortuosity [53] 3D vuggy limestone 0.16-0.81Permeability [21] 3D samples of volcanic tuff, glass bead column, sandpack, sandstones, carbonates 3D samples generated by level-set percolation Real: 0.43, 0.6 Generated: 0.35, 0.6 Permeability, hydraulic tortuosity [12] 3D carbonate, sandpacks, sandstone 0.17-0.38Permeability [22] 3D samples generated by spherical/nonspherical grain models 0.05-0.55 Permeability, hydraulic tortuosity [13] 2D randomly distributed freely overlapping squares 0.367-0.99Hydraulic tortuosity 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.

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 u by where k is the permeability tensor.The other quantities are fluid viscosity and density  and , macroscale pressure , and gravitational acceleration g.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.
(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].

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 is where   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  2 ℎ 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.

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.

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.Isotropic system, 0 <  < 0.5  =  −0.4 [31] Fine-grained uncemented sediments  = (1 − ln ( 2 )) [56] Random, partial overlapping shapes  = 1 1 −  (1 − )  = shape factor [46] Voronoi channel geometries, 0.2 <  < 0.6  = (0.75 −1.0819 ) 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 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 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.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.

Methods to Compute Effective Properties
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  = 1, 2, 3 and  = 1, 2, 3 and repeated indices imply summation,  is pressure,  is fluid viscosity, V  is the pore space fluid velocity,  is fluid density, and   is gravitational acceleration.The boundary condition at the fluid-solid interface is given by which 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 V , 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 0 ≤  ≤ 1.A multiscale expansion is used to write   and V , in terms of their leading and higher order terms as follows: V , = V  (,   ) = V (0)  (, ) + V (1)   (, ) +  2 V (2)  (, ) + ⋅ ⋅ ⋅ .
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 are such that   () and   () satisfy with the fluid-solid boundary condition given by In 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 (  = 1 for  = , and   = 0 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 V (2) (, ) in the cell is found to be divergence free; ∇ ⋅ ⟨V (2) (, )⟩  = 0.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 as where  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 V = (V  , V  )  and the pressure  such that where 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 ( +1/2 ,  +1/2 ) is ( +1/2 ,  +1/2 ) approximated by  +1/2,+1/2 .Similarly, we approximate the -velocity V  at each -edge (fluid) center V  (  ,  +1/2 ) by V  ,+1/2 , and the -velocity 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 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 computergenerated media samples are 2-dimensional (2D); thus the set is solved with different forcings   where  = 1, 2 and  = 1, 2. 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 k and k * for the permeability tensors in coordinate systems {x} and {x * }, 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 is where ⟨⋅⟩ Ω is the spatial average in the unit cell domain Ω, V mag is the magnitude of the fluid velocities (i.e., the fluid speed), and V  is the velocity component in the  direction.
The spatial average of V  in Ω is where  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 V mag 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  2 ℎ 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 V  and V  explicitly.Instead, during the permeability computation for our 2D pore geometries, we obtain the rescaled pore space velocity fields   in the cell, where  = 1, 2 and  = 1, 2. It can be seen that  11 and  21 are related to V  and V  , respectively, for flow that is driven by g = (ĝ, 0) and that  12 and  22 are related to V  and V  , respectively, for flow that is driven by g = (0, ĝ).Thus, it is mathematically equivalent to compute the hydraulic tortuosities by where Ω 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   = 0 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 +/2 fl θ.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  is where 2 and where   1 ,   2 , and   are obtained by vector transformations from   . ℎ in the other principal direction θ is similarly obtained.The notation we use in our results section is  ℎ,max =  ℎ, and  ℎ, =  ℎ, θ.

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 by where V , 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 expansion and ( 8) are substituted into the governing transport equations (30).Upon collecting like order terms, the second-order term of  is where  (0) () is the macroscale concentration and   (, ) where (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  eff *  is given by Intuitively, the more tortuous the porous media, the slower the rate of effective diffusion; thus a commonly used relationship between effective and molecular diffusivity is where  , 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][43][44][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 by where   =  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  ,max and  , corresponding to the principal directions.The rotation required to obtain the principal components of k * 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).

Results
While real geological formations are typically of low porosity ( < 0.3), we have generated several idealized geometries within the porosity range of 0.33 <  < 1.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 tortuosityporosity 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.

Permeability and Hydraulic
Tortuosity.The full permeability tensor was computed for 0.33 <  < 1, 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 Table 4: Input parameters used to generate different pore-geometries (or different -cells) and their corresponding porosity ranges.In this table, cell refers to a cell in the mesh used to discretize the -cell.

Geometry
Initial mesh dimensions Inputs Porosity range In-line array of shapes (see Figure 1 our permeability data on increasingly finer meshes until the permeability converged within a 1% relative error.To reach this convergence criteria after starting with a mesh size of 50 × 50, the lowest porosity cell required a refined mesh of 500 × 500, and the highest porosity cell required 300 × 300.
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 0.33 <  < 1, with the same convergence criteria used to obtain permeability.Figure 3 illustrates the fluid speed and direction for a unit cell of  = 0.71, and the hydraulic tortuosity in both  and  directions were computed to be  ℎ, =  ℎ, = 1.0185.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.

Diffusive Tortuosity.
We obtained the full diffusive tortuosity tensor on increasingly refined meshes until the tortuosity converged within a 0.1% relative error.Starting with a mesh of 50 × 50, the lowest porosity required a refined mesh of 700 × 700, while the highest porosity required 250 × 250 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.described in Section 3.3.The resulting diffusive tortuosity tensor is where   = (   ) −1 .Notice the tensor is symmetric and diagonal (i.e.,  , =  , = 0 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 ∇ 2  = 0 with Dirichlet boundary conditions on the right and left boundaries such that the concentration gradient across the unit cell is one.

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  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 (0,  , ) is tabulated in Table 5.These lengths are computed by   =   /, where   is the length of the th streamline and  = 1 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  = 2 −  (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 (  = 1 − /4 ≈ 0.21) 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.

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, / = 1 and   /  = 1.By varying the solid square length, a range of porosities are obtained.

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 100×50, the meshes required refinement up to a maximum of 700×350 in order to meet the convergence criteria of  rel error < 1%.
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  = 0.73 are shown in Figure 11 and the resulting hydraulic tortuosities were  ℎ, = 1.3539 and  ℎ, = 1.0567.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   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).

Diffusive Tortuosity.
We obtained the full diffusive tortuosity tensor for this staggered geometry using a convergence criteria of  ,rel error < 0.1%.Starting with a uniform mesh of 100 × 50, the lowest porosity required a refined mesh of 1000 × 500 while the highest porosity required 200 × 100 in order to meet the convergence criteria.Results for a   staggered-array unit cell of  = 0.73 are shown in Figure 12, and the diffusive tortuosity tensor is where the tensor is symmetric and diagonal (i.e.,  , =  , = 0 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.

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 / = 1,   /  = 1), 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 D eff / rather than   explicitly, it may not have been easily apparent to them that one of their values corresponded to   < 1. (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 ( < 0.6).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  ℎ, ≈ 1 for all porosities while  ℎ, increases towards the lower porosities.

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 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 anisotropic-related statistical errors.This same constraint was also respected in Duda et al. [13].Our geometries ranged from (effective) porosities of 0.45 <  < 0.99.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.

Permeability and Hydraulic
Tortuosity.We computed permeability using a convergence criteria of  rel error < 1%.
Starting with an initial mesh of 100 × 100, the lowest porosity required a refined mesh of 400 × 400 and the highest porosity required 200 × 200 to meet the convergence criteria.The full permeability tensors (in coordinate system {x}) for these pore geometries were symmetric and nondiagonal; thus a transformation was done to obtain the diagonal tensor in the principal coordinate system {x * }.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).
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  = 0.7, results for hydraulic flow are shown in Figure 17 and the computed hydraulic tortuosities are  ℎ, = 1.3727 and  ℎ, = 1.3393.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 V mag > 5 × 10 −5 ) exist in the geometry, while the rest of the pore space is characterized by moderate (1 × 10 −5 < V mag < 5 × 10 −5 ) to slow (V mag < 1 × 10 −5 ) moving fluid.The high speeds  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.

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 0.7 <  < 0.90, which was similarly reported in Duda et al. [13].

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  rel error < 3%, 2%, and 1% were used for porosity ranges of 0.45 <  < 0.58, 0.58 <  < 0.7, and 0.7 <  < 0.99, respectively.The geometries were generated using a mesh of 100 × 100, and the convergence criteria required a refinement of 2000 × 2000 for the lowest porosity and 200 × 200 for the highest porosity.Results for a geometry of  = 0.7 are shown in Figure 21, and the resulting diffusive tortuosity tensor is Notice 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.

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  < 0.6 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.
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, k * 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).)

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 k; 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 k 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   Geofluids 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, k 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.

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. where   =   .Note that diffusive tortuosity is the inverse of the diffusive tortuosity factor, that is,   = (   ) −1 .Recalling (37), the diffusive tortuosity factor is related to the effective diffusion coefficient by where D eff = D eff * .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.

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

Figure 1 :
Figure 1: Macroscale (a) and microscale (b) domains used in the homogenization problem.The solid structure in the -cell illustrates the geometry of our first example: in-line array of solid shapes (i.e., circles or squares).

) 50 ×Figure 3 :
Figure 3: Fluid flow past a solid circle,  = 0.71, driven by   .Nondimensional fluid speed (i.e., magnitude of velocities) is indicated by color, and flow path is indicated by flow lines (direction is left to right).Speed profile and flow path driven by   are same as shown but rotated 90 degrees.

Figure 4 :
Figure 4: Number of grid cells required to reach converged   (using a convergence criteria of 0.1% relative error) for select porosities of in-line array of circles.

Figure 5 (Figure 5 :Figure 6 :
Figure 5: Fields for a unit cell of  = 0.71, comprised of a solid circle: (a) solution to the homogenized problem, (b) concentration field transformed from   , and (c) magnitude of concentration gradient √(∇    ) 2 + (∇    ) 2 (diffusion speed) and direction.The fields with respect to  are the same as shown above but rotated by 90 degrees.
(c)   appear to indicate that the flow fields are relatively analogous.However, the fluid speeds shown in Figure3are such that the fastest moving flow exists along the nontortuous pathway, while the fastest passive diffusive transport (or highest concentration gradient) shown in Figure5(c) exists along

Figure 10 :
Figure10: Anisotropic ratios in staggered-array of squares.For all porosities, the anisotropic ratios for hydraulic and diffusive tortuosity are less than one, indicating that the hydraulic flow and passive diffusive transport are more tortuous in the  direction than in the  direction.The hydraulic flow is more permeable in the  direction.
Flow in  direction

Figure 11 :
Figure 11: Fluid flow through a staggered-array of squares,  = 0.73.Fluid speed (i.e., magnitude of velocities) is indicated by color, and direction is indicated by flow lines.Fluid speed is nondimensional.

Figure 12 :
Figure 12: Fields for a staggered-array unit cell of  = 0.73: (a, d) solution to the homogenization problem, (b, e) concentration field transformed from   (assuming concentration gradient across unit cell is one), and (c, f) magnitude of concentration gradient (diffusion speed) and direction.

Figure 13 :
Figure13: Tortuosity-porosity trends for staggered-array of squares.Simulated data is plotted against (diffusive) data from literature that employed finite difference (FD) and boundary element (BE) methods.In this geometry, hydraulic tortuosity in the  direction is a weak function of porosity.

Figure 14 :
Figure 14: Various pore structures for randomly distributed squares geometry.Any fluid area isolated from flow path is converted to solid; thus  is effective porosity of geometry.
Figure 17: Flow (driven by   ) through randomly distributed squares,  = 0.70.Nondimensional fluid speed is indicated by color.Zero speed at solid sites is kept as dark blue for clarity of illustration (while these fluid speeds do not actually exist).

Figure 18 :
Figure 18: Hydraulic tortuosity-porosity trend for randomly distributed squares.The flow pathways are more tortuous as geometry is comprised of more solid squares (i.e., lower porosity).

Figure 19 :
Figure19: Specific surface area for randomly distributed squares, over the range of porosities studied.Specific surface area increases linearly as porosity decreases; however this trend deviates at lower porosities due to the large nonuniform shapes introduced into the lower porosity geometries.
Porosityfit for k max and  h,max data fit for k t and  h,t data

Figure 20 :
Figure20: Shape factor,   , used to fit simulated permeability and hydraulic tortuosity data to Kozeny-Carman equation for geometries comprised of randomly distributed squares.  appears to be constant between 0.7 <  < 0.9.

Figure 21 :
Figure 21: Fields for geometry of randomly distributed squares,  = 0.7: (a) solution to the homogenization problem, (b) concentration field transformed from   (assuming concentration gradient across unit cell is one), and (c, d) magnitude of concentration gradient (diffusion speed).

Figure 22 :
Figure 22: Diffusive tortuosity-porosity trend for randomly distributed squares.The full diffusive tortuosity tensor contains nonzero offdiagonal components and is diagonalized to obtain the principal diffusive tortuosity components,  ,max and  , .{x} and {x * } are related by rotation angle   , which is plotted in Figure 16(a).The anisotropic ratio of the principal components is plotted in Figure 16(b).

Figure 24 :Figure 25 :
Figure 24: Induced anisotropy in lower porosity geometries comprised of randomly distributed squares.

Table 1 :
Recent work using real or computer-generated porous media samples: comparison of sample types, porosities, and computed properties.

Table 2 :
Past work on hydraulic tortuosity.

Table 3 :
Past work on diffusive tortuosity.

Table 5 :
Hydraulic and diffusive streamline comparison: in-line array of circles,  = 0.71.Flow is from left to right.(0,  , ) indicate the starting location of the th streamline.
the circle's rounded fluid-solid boundary.This difference in speeds (or streamline weights) leads to a quantitative difference between these two tortuosity types.

Table 6 :
Data comparison for diffusive tortuosity. , 3, and for a 2D problem the tensor components are computed by [  ,  ,  ,  , ]