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

Authors: Rebecca Allen; Shuyu Sun
Journal: Geofluids
 (2017)

---

## Abstract

We compute effective properties (i.e., permeability, hydraulic tortuosity, and diffusive tortuosity) of three different digital porous media samples, including in-line array of uniform shapes, staggered-array of squares, and randomly distributed squares. The permeability and hydraulic tortuosity are computed by solving a set of rescaled Stokes equations obtained by homogenization, and the diffusive tortuosity is computed by solving a homogenization problem given for the effective diffusion coefficient that is inversely related to diffusive tortuosity. We find that hydraulic and diffusive tortuosity can be quantitatively different by up to a factor of ten in the same pore geometry, which indicates that these tortuosity terms cannot be used interchangeably. We also find that when a pore geometry is characterized by an anisotropic permeability, the diffusive tortuosity (and correspondingly the effective diffusion coefficient) can also be anisotropic. This finding has important implications for reservoir-scale modeling of flow and transport, as it is more realistic to account for the anisotropy ofboth the permeability and the effective diffusion coefficient.


---

## Body

## 1. Introduction

When modeling subsurface flow and transport in a reservoir, it is common practice to represent the geological formation with effective properties instead of resolving the precise location of fluid and solid phase. Effective properties such as permeability, diffusivity, and tortuosity can be computed from pore-scale modeling of fluid flow and transport through a porous media sample (i.e., a core sample) taken from a geological formation. Permeability and effective diffusivity are the well-known and important properties used in Darcy- and reservoir-scale flow and transport equations, while tortuosity is given less importance and subject to multiple definitions. Tortuosity is generally defined as a ratio between the effective path traveled by a species and the unit length of the domain. As noted in [1], different types of tortuosity (geometric, hydraulic, diffusive, dispersive, and electrical) are distinguished by the type of species transport under consideration. Hydraulic tortuosity refers to the effective path traveled by a fluid particle driven by a force, while diffusive tortuosity refers to the effective path traveled by a species driven by molecular diffusion [1].Pore-scale modeling is performed on digital samples of porous media, which can be obtained using imaging and conversion techniques as those described in [2, 3], or by using an algorithm that generates samples which possess the same statistical characteristics of a real porous media sample [4–6]. Table 1 presents a summary of studies which used eitherreal orcomputer-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 ofrecent works on the topic, and that pore-scale simulation to obtain effective properties dates back to Cancelliere et al. [7] and possibly earlier.Table 1
Recent work using real or computer-generated porous media samples: comparison of sample types, porosities, and computed properties.
Study Real sample Computer-generated sample Porosity range Effective properties computed [16]  3D fibrous material (membrane) 0.45–0.9 permeability, 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.81 Permeability  [21] 3D samples of volcanic tuff, glass bead column, sandpack, sandstones, carbonates 3D samples generated by level-set percolation Real: 0.43, 0.6Generated: 0.35, 0.6 Permeability, hydraulic tortuosity  [12] 3D carbonate, sandpacks, sandstone  0.17–0.38 Permeability  [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.99 Hydraulic tortuosityIn regard to obtaining effective properties at the pore-scale, we are motivated to distinguish between two different types of tortuosity (namely, hydraulic and diffusive) that have been identified. While a few studies have compared tortuosity definitions and have shown they are indeed quantitatively different within the same pore structure [8, 9], other studies continue to use these definitions interchangeably or do not distinguish the type of tortuosity they are computing. The usefulness of hydraulic tortuosity becomes apparent when formulating a prediction for permeability, such as the Kozeny-Carman equation. Also, tortuosity can be seen as a more intuitive quantity for fluid flow and transport through a pore geometry than permeability or the effective diffusion coefficient.The outline of this work is as follows: first we will conduct a literature review on the methods commonly used to compute the effective properties of digital porous media; then we will present the details of these methods we employed herein, followed by three example geometries (in-line array of uniform shapes, staggered-array of squares, and randomly distributed squares) and our obtained trends between porosity and the effective properties. We will conclude with the main findings of this work.


## 2. Literature Review on Effective Properties

### 2.1. Permeability

From the list given in Table1, 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(1)u=-kμ·∇P+ρg,where k is the permeability tensor. The other quantities are fluid viscosity and density μ and ρ, macroscale pressure P, 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 Table1, the general approach to compute the permeability field of the porous media sample is as follows:(1)
Obtain pore-scale geometries: converted from real media sample images or generated by an algorithm that may or may not use media parameters such as grain size distribution and so on.(2)
Solve Stokes flow in pore space of media: different solvers have been used (finite difference [12], finite element, lattice Boltzmann [13–16], explicit jump method [15], etc.), from academically developed codes, commercial software, and open-source software.(3)
Compute permeability: assuming Darcy’s equation is valid for the pore structure and flow field results, permeability is computed based on volume averages of velocity and pressure field.A further step is to fit the permeability data to a formula which proposes a relationship between permeability and other porous media properties, known as the Kozeny-Carman (KC) equation [17, 18].


### 2.2. Kozeny-Carman Equation

An idealized pore geometry of parallel cylindrical channels was used to formulate the Kozeny-Carman (KC) equation [17]. This equation proposes a relationship between porosity ϕ and permeability k, as well as other parameters such as hydraulic tortuosity τh and specific surface area S, and is(2)k=ϕ3βkτh2S2,where βk 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 τh in (2) is the ratio between effective path traveled by the fluid particle and the length of the sample: τh=Le/Ls. 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 τh2 as the hydraulic tortuosity as similarly pointed out in [19]. The specific surface area is the ratio of fluid-solid interface area Af-s to the total volume in the domain:(3)S=Aporesf-sVtot.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 k=kii and τh=τh,ii 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 βk 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] measureddiffusive 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 computedhydraulic tortuosity. Also, Sun et al. [27] applied the method of homogenization to the diffusive transport equation in order to compute the effective diffusion coefficient of periodic unit cells and the corresponding tortuosity. What they did not point out was that their approach gave them the diffusive tortuosity, while some of the trends they compared their results to were with respect to hydraulic tortuosity. Unless tortuosity types are quantitatively identical in the same class of porous media, use of a trend that is specific to one type of tortuosity should not immediately be used to describe the trend of another type of tortuosity.A few studies [8, 9, 28] have focused on quantifying the difference between two specific tortuosities in the same structure. For example, hydraulic and electrical tortuosity were compared in both David [28] and Zhang and Knackstedt [9]. David [28] used a network as an analog for the media pore space, and the study focused on quantifying the tortuosity of the critical or preferential pathway, not the tortuosity of the entire flow field. Zhang and Knackstedt [9] computed tortuosity using the entire flow field which they obtained by numerical simulations for fluid flow (via lattice Gas Automata) and electrical current (via finite difference). Both of these works found that the hydraulic tortuosity was higher than electrical tortuosity in the same structure (up to an order of magnitude in low porosity configurations [9]). Once again, it was emphasized in Ghanbarian et al.’s [8] review that due to the difference between tortuosity types their models are not interchangeable.Regarding diffusive tortuosity, both Quintard et al. [29] and Valdés-Parada et al. [1] stated that its quantification becomes more complex when more than just passive diffusion is occurring in the system. Valdés-Parada et al. [1] showed how the quantification of diffusive tortuosity can be different depending on the consideration of different mass transport terms (i.e., passive diffusion only, or diffusion with convection, reaction, or hydrodynamic dispersion). They concluded that the consideration of passive diffusion leads to the only appropriate definition of tortuosity. Their work focused on mass transport of species through fluid and did not measure hydraulic tortuosity.


#### 2.3.2. Various Tortuosity-Porosity Trends

Many authors have theoretically or empirically derived tortuosity as a function of porosity (see Tables2 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.Table 2
Past work on hydraulic tortuosity.
Study Samples considered Tortuosity–porosity fit Notes [54] Packed beds τ = 1 + p ln ⁡ 1 ϕ p is function of particle shape [19] 2D random overlapping monosized squares,0.5<ϕ<1 τ = 1 + 0.8 1 - ϕ [26] 2D random overlapping monosized squares,0.4<ϕ<1 τ = 1 + 0.65 1 - ϕ ϕ - ϕ c 0.19 ϕ c is critical porosity [24] 2D random overlapping monosized squares τ ∝ 1 + R S ϕ R is hydraulic radius, S is specific surface area [13] 2D freely overlapping squares τ ∝ 1 + 1 - ϕ 1 / 2Table 3
Past work on diffusive tortuosity.
Study Samples considered Tortuosity–porosity fit Notes [48]∗ Array of cylinders (2D) τ = 2 - ϕ [50]∗∗ Diffusion of electrolytes in membrane τ = 2 - ϕ ϕ 2 [52] Bed of uniform spheres (applicable to overlapping, nonuniform spheres) τ = 1 - 1 2 ln ⁡ ϕ [43] Isotropic system,0<ϕ<0.5 τ = ϕ - 0.4 [31] Fine-grained uncemented sediments τ = 1 - ln ⁡ ϕ 2 1 / 2 [55] Isotropic representative unit cell τ = ϕ 1 - 1 - ϕ m 2D:m=12, 3D: m=23 [56] Random, partial overlapping shapes τ = 1 1 - α 1 - ϕ α = shape factor [46] Voronoi channel geometries,0.2<ϕ<0.6 τ = 0.75 ϕ - 1.0819 1 / 2 ∗ As referenced in [57]. ∗ ∗ As referenced in [30, 31].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 thatany definition of tortuosity should not be considered as a function of porosity but rather a function of the pore geometry only. The idea behind this statement seems to rest in the fact that there is no universally agreed upon trend between tortuosity and porosity as pointed out by Matyka et al. [24] and that the geometrical features of a pore structure should have more influence on the overall tortuosity of flow or transport through the geometry than the value of porosity. However, since it is possible that a unique relationship between tortuosity and porosity can indeed exist for special classes of porous media [24], the development of tortuosity-porosity trends remains to be an insightful topic of research.Also, Liu and Kitanidis [35] stated that tortuosity is a tensorial property. Despite this statement, the majority of work does not report directional-specific quantities which may leave the reader with the impression that tortuosity is a scalar quantity. A scalar quantity could be appropriate in pore geometries that exhibit an anisotropic microstructure while still exhibiting an isotropic macrostructure. For example, pore geometry composed of randomly distributed squares is geometrically anisotropic at the microscale; however the effective properties computed for the REV can be isotropic in nature. In this work, we demonstrate this point by example.


## 2.1. Permeability

From the list given in Table1, 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(1)u=-kμ·∇P+ρg,where k is the permeability tensor. The other quantities are fluid viscosity and density μ and ρ, macroscale pressure P, 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 Table1, the general approach to compute the permeability field of the porous media sample is as follows:(1)
Obtain pore-scale geometries: converted from real media sample images or generated by an algorithm that may or may not use media parameters such as grain size distribution and so on.(2)
Solve Stokes flow in pore space of media: different solvers have been used (finite difference [12], finite element, lattice Boltzmann [13–16], explicit jump method [15], etc.), from academically developed codes, commercial software, and open-source software.(3)
Compute permeability: assuming Darcy’s equation is valid for the pore structure and flow field results, permeability is computed based on volume averages of velocity and pressure field.A further step is to fit the permeability data to a formula which proposes a relationship between permeability and other porous media properties, known as the Kozeny-Carman (KC) equation [17, 18].


## 2.2. Kozeny-Carman Equation

An idealized pore geometry of parallel cylindrical channels was used to formulate the Kozeny-Carman (KC) equation [17]. This equation proposes a relationship between porosity ϕ and permeability k, as well as other parameters such as hydraulic tortuosity τh and specific surface area S, and is(2)k=ϕ3βkτh2S2,where βk 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 τh in (2) is the ratio between effective path traveled by the fluid particle and the length of the sample: τh=Le/Ls. 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 τh2 as the hydraulic tortuosity as similarly pointed out in [19]. The specific surface area is the ratio of fluid-solid interface area Af-s to the total volume in the domain:(3)S=Aporesf-sVtot.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 k=kii and τh=τh,ii 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 βk 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] measureddiffusive 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 computedhydraulic tortuosity. Also, Sun et al. [27] applied the method of homogenization to the diffusive transport equation in order to compute the effective diffusion coefficient of periodic unit cells and the corresponding tortuosity. What they did not point out was that their approach gave them the diffusive tortuosity, while some of the trends they compared their results to were with respect to hydraulic tortuosity. Unless tortuosity types are quantitatively identical in the same class of porous media, use of a trend that is specific to one type of tortuosity should not immediately be used to describe the trend of another type of tortuosity.A few studies [8, 9, 28] have focused on quantifying the difference between two specific tortuosities in the same structure. For example, hydraulic and electrical tortuosity were compared in both David [28] and Zhang and Knackstedt [9]. David [28] used a network as an analog for the media pore space, and the study focused on quantifying the tortuosity of the critical or preferential pathway, not the tortuosity of the entire flow field. Zhang and Knackstedt [9] computed tortuosity using the entire flow field which they obtained by numerical simulations for fluid flow (via lattice Gas Automata) and electrical current (via finite difference). Both of these works found that the hydraulic tortuosity was higher than electrical tortuosity in the same structure (up to an order of magnitude in low porosity configurations [9]). Once again, it was emphasized in Ghanbarian et al.’s [8] review that due to the difference between tortuosity types their models are not interchangeable.Regarding diffusive tortuosity, both Quintard et al. [29] and Valdés-Parada et al. [1] stated that its quantification becomes more complex when more than just passive diffusion is occurring in the system. Valdés-Parada et al. [1] showed how the quantification of diffusive tortuosity can be different depending on the consideration of different mass transport terms (i.e., passive diffusion only, or diffusion with convection, reaction, or hydrodynamic dispersion). They concluded that the consideration of passive diffusion leads to the only appropriate definition of tortuosity. Their work focused on mass transport of species through fluid and did not measure hydraulic tortuosity.


### 2.3.2. Various Tortuosity-Porosity Trends

Many authors have theoretically or empirically derived tortuosity as a function of porosity (see Tables2 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.Table 2
Past work on hydraulic tortuosity.
Study Samples considered Tortuosity–porosity fit Notes [54] Packed beds τ = 1 + p ln ⁡ 1 ϕ p is function of particle shape [19] 2D random overlapping monosized squares,0.5<ϕ<1 τ = 1 + 0.8 1 - ϕ [26] 2D random overlapping monosized squares,0.4<ϕ<1 τ = 1 + 0.65 1 - ϕ ϕ - ϕ c 0.19 ϕ c is critical porosity [24] 2D random overlapping monosized squares τ ∝ 1 + R S ϕ R is hydraulic radius, S is specific surface area [13] 2D freely overlapping squares τ ∝ 1 + 1 - ϕ 1 / 2Table 3
Past work on diffusive tortuosity.
Study Samples considered Tortuosity–porosity fit Notes [48]∗ Array of cylinders (2D) τ = 2 - ϕ [50]∗∗ Diffusion of electrolytes in membrane τ = 2 - ϕ ϕ 2 [52] Bed of uniform spheres (applicable to overlapping, nonuniform spheres) τ = 1 - 1 2 ln ⁡ ϕ [43] Isotropic system,0<ϕ<0.5 τ = ϕ - 0.4 [31] Fine-grained uncemented sediments τ = 1 - ln ⁡ ϕ 2 1 / 2 [55] Isotropic representative unit cell τ = ϕ 1 - 1 - ϕ m 2D:m=12, 3D: m=23 [56] Random, partial overlapping shapes τ = 1 1 - α 1 - ϕ α = shape factor [46] Voronoi channel geometries,0.2<ϕ<0.6 τ = 0.75 ϕ - 1.0819 1 / 2 ∗ As referenced in [57]. ∗ ∗ As referenced in [30, 31].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 thatany definition of tortuosity should not be considered as a function of porosity but rather a function of the pore geometry only. The idea behind this statement seems to rest in the fact that there is no universally agreed upon trend between tortuosity and porosity as pointed out by Matyka et al. [24] and that the geometrical features of a pore structure should have more influence on the overall tortuosity of flow or transport through the geometry than the value of porosity. However, since it is possible that a unique relationship between tortuosity and porosity can indeed exist for special classes of porous media [24], the development of tortuosity-porosity trends remains to be an insightful topic of research.Also, Liu and Kitanidis [35] stated that tortuosity is a tensorial property. Despite this statement, the majority of work does not report directional-specific quantities which may leave the reader with the impression that tortuosity is a scalar quantity. A scalar quantity could be appropriate in pore geometries that exhibit an anisotropic microstructure while still exhibiting an isotropic macrostructure. For example, pore geometry composed of randomly distributed squares is geometrically anisotropic at the microscale; however the effective properties computed for the REV can be isotropic in nature. In this work, we demonstrate this point by example.


## 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] measureddiffusive 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 computedhydraulic tortuosity. Also, Sun et al. [27] applied the method of homogenization to the diffusive transport equation in order to compute the effective diffusion coefficient of periodic unit cells and the corresponding tortuosity. What they did not point out was that their approach gave them the diffusive tortuosity, while some of the trends they compared their results to were with respect to hydraulic tortuosity. Unless tortuosity types are quantitatively identical in the same class of porous media, use of a trend that is specific to one type of tortuosity should not immediately be used to describe the trend of another type of tortuosity.A few studies [8, 9, 28] have focused on quantifying the difference between two specific tortuosities in the same structure. For example, hydraulic and electrical tortuosity were compared in both David [28] and Zhang and Knackstedt [9]. David [28] used a network as an analog for the media pore space, and the study focused on quantifying the tortuosity of the critical or preferential pathway, not the tortuosity of the entire flow field. Zhang and Knackstedt [9] computed tortuosity using the entire flow field which they obtained by numerical simulations for fluid flow (via lattice Gas Automata) and electrical current (via finite difference). Both of these works found that the hydraulic tortuosity was higher than electrical tortuosity in the same structure (up to an order of magnitude in low porosity configurations [9]). Once again, it was emphasized in Ghanbarian et al.’s [8] review that due to the difference between tortuosity types their models are not interchangeable.Regarding diffusive tortuosity, both Quintard et al. [29] and Valdés-Parada et al. [1] stated that its quantification becomes more complex when more than just passive diffusion is occurring in the system. Valdés-Parada et al. [1] showed how the quantification of diffusive tortuosity can be different depending on the consideration of different mass transport terms (i.e., passive diffusion only, or diffusion with convection, reaction, or hydrodynamic dispersion). They concluded that the consideration of passive diffusion leads to the only appropriate definition of tortuosity. Their work focused on mass transport of species through fluid and did not measure hydraulic tortuosity.


## 2.3.2. Various Tortuosity-Porosity Trends

Many authors have theoretically or empirically derived tortuosity as a function of porosity (see Tables2 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.Table 2
Past work on hydraulic tortuosity.
Study Samples considered Tortuosity–porosity fit Notes [54] Packed beds τ = 1 + p ln ⁡ 1 ϕ p is function of particle shape [19] 2D random overlapping monosized squares,0.5<ϕ<1 τ = 1 + 0.8 1 - ϕ [26] 2D random overlapping monosized squares,0.4<ϕ<1 τ = 1 + 0.65 1 - ϕ ϕ - ϕ c 0.19 ϕ c is critical porosity [24] 2D random overlapping monosized squares τ ∝ 1 + R S ϕ R is hydraulic radius, S is specific surface area [13] 2D freely overlapping squares τ ∝ 1 + 1 - ϕ 1 / 2Table 3
Past work on diffusive tortuosity.
Study Samples considered Tortuosity–porosity fit Notes [48]∗ Array of cylinders (2D) τ = 2 - ϕ [50]∗∗ Diffusion of electrolytes in membrane τ = 2 - ϕ ϕ 2 [52] Bed of uniform spheres (applicable to overlapping, nonuniform spheres) τ = 1 - 1 2 ln ⁡ ϕ [43] Isotropic system,0<ϕ<0.5 τ = ϕ - 0.4 [31] Fine-grained uncemented sediments τ = 1 - ln ⁡ ϕ 2 1 / 2 [55] Isotropic representative unit cell τ = ϕ 1 - 1 - ϕ m 2D:m=12, 3D: m=23 [56] Random, partial overlapping shapes τ = 1 1 - α 1 - ϕ α = shape factor [46] Voronoi channel geometries,0.2<ϕ<0.6 τ = 0.75 ϕ - 1.0819 1 / 2 ∗ As referenced in [57]. ∗ ∗ As referenced in [30, 31].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 thatany definition of tortuosity should not be considered as a function of porosity but rather a function of the pore geometry only. The idea behind this statement seems to rest in the fact that there is no universally agreed upon trend between tortuosity and porosity as pointed out by Matyka et al. [24] and that the geometrical features of a pore structure should have more influence on the overall tortuosity of flow or transport through the geometry than the value of porosity. However, since it is possible that a unique relationship between tortuosity and porosity can indeed exist for special classes of porous media [24], the development of tortuosity-porosity trends remains to be an insightful topic of research.Also, Liu and Kitanidis [35] stated that tortuosity is a tensorial property. Despite this statement, the majority of work does not report directional-specific quantities which may leave the reader with the impression that tortuosity is a scalar quantity. A scalar quantity could be appropriate in pore geometries that exhibit an anisotropic microstructure while still exhibiting an isotropic macrostructure. For example, pore geometry composed of randomly distributed squares is geometrically anisotropic at the microscale; however the effective properties computed for the REV can be isotropic in nature. In this work, we demonstrate this point by example.


## 3. Methods to Compute Effective Properties

### 3.1. Permeability

The method of homogenization has been used to mathematically derive Darcy’s equation from the Stokes equations. Through this derivation, the macroscopic property known as permeabilitykij is defined as the spatial average of arescaled pore space velocity Nij(y). 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 Figure1, 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:(4)-∂∂xipε+μ∂∂xj∂∂xjvε,i=-ρgj,(5)∂∂xivε,i=0,where i=1,2,3 and j=1,2,3 and repeated indices imply summation, p is pressure, μ is fluid viscosity, vi is the pore space fluid velocity, ρ is fluid density, and gj is gravitational acceleration. The boundary condition at the fluid-solid interface is given by(6)vi=0,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].Figure 1
Macroscale (a) and microscale (b) domains used in the homogenization problem. The solid structure in theY-cell illustrates the geometry of our first example: in-line array of solid shapes (i.e., circles or squares).In the Stokes equations,pε and vε,i can vary within the pore space and thus are functions of the macroscale (x) and microscale (ε); however μ, g, and ρ are treated as constants and are thus functions of the macroscale (x) only. By defining a fast variable y=x/ε, the cell of unit length ε is scaled into the Y-cell (shown in Figure 1(b)), where 0≤y≤1. A multiscale expansion is used to write pε and vε,i in terms of their leading and higher order terms as follows:(7)pε=px,xε=p0x,y+εp1x,y+⋯,(8)vε,i=vix,xε=vi0x,y+εvi1x,y+ε2vi2x,y+⋯.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(9)p1x,y=Qjy-ρgj+∂∂xjp0x,v2x,y=Nijyμ-ρgj+∂∂xjp0x,such that Qj(y) and Nij(y) satisfy(10)-∂∂yiQjy+∂∂yk∂∂ykNijy=δij,(11)∂∂yiNijy=0,with the fluid-solid boundary condition given by(12)Nijy=0.In the above formulation, Qj(y) and Nij(y) are functions of the microscale only and are independent of macroscale quantities μ, ρ, and gj. 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 Y-cell. In (10), δij is the Kronecker delta function (δij=1 for i=j, and δij=0 for i≠j).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 ofv(2)(x,y) in the cell is found to be divergence free; ∇·v2x,yy=0. This implies that the macroscopic flow equation is satisfied by Darcy’s equation (1), that is,(13)∂∂xikijμ∂∂xjpx-ρgj=0⟶∂∂xiuix=0,where p(x) is the macroscale pressure and ui(x) is the macroscale velocity, when the full permeability tensor for an n-dimensional problem is defined as(14)kij≔Nijyy=1Y∫YNijydY,where Y is the volume of the Y-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 Y-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=(vx,vy)T and the pressure p such that(15)-μΔv+∇p=f,inΩ≔0,Lx×0,Ly,(16)∇·v=0,inΩ,(17)v=0,on⁡∂Ω,where the velocity has a unique solution while the pressure is unique up to an additive constant. Here Lx and Ly denote the length of the rectangular domain Ω in the x- and y-directions, respectively.To construct a staggered grid finite difference scheme for the Stokes problem, we introduce a possibly nonuniform grid ofΩ as follows: (18)0=x0<x1<⋯<xM=Lx,0=y0<y1<⋯<yN=Ly.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 (xi+1/2,yj+1/2) is p(xi+1/2,yj+1/2) approximated by pi+1/2,j+1/2. Similarly, we approximate the x-velocity vx at each x-edge (fluid) center vx(xi,yj+1/2) by vi,j+1/2x, and the y-velocity vy at each y-edge (fluid) center vy(xi-1/2,yj) by vi-1/2,jy.By integrating thex-component of (15) over cell Ωi,j+1/2=(xi-1/2,xi+1/2)×(yj,yj+1), we have(19)-∫yjyj+1∂vx∂xxi+1/2,y-∂vx∂xxi-1/2,ydy+∫xi-1/2xi+1/2∂vx∂yx,yj+1-∂vx∂yx,yjdx+∫yjyj+1pxi+1/2,y-pxi-1/2,ydy=∫Ωi,j+1/2fxx,ydxdy.We then approximate those edge integrals by midpoint values and the volume integral by its center value and divide both sides byΩi,j+1/2=(xi+1/2-xi-1/2)(yj+1-yj).By introducing the following forward and backward difference operators:(20)Δx+vi,j+1/2x=vi+1,j+1/2x-vi,j+1/2xxi+1-xi,Δx-pi+1/2,j+1/2=pi+1/2,j+1/2-pi-1/2,j+1/2xi+1/2-xi-1/2and similar notations for operations in the y-direction, we can express the results as(21)-Δx+vi,j+1/2x-Δx-vi,j+1/2xxi+1/2-xi-1/2-Δy+vi,j+1/2x-Δy-vi,j+1/2xyj+1-yj+Δx-pi+1/2,j+1/2=fi,j+1/2x.Similarly, by integrating they-component of (15) over cell Ωi-1/2,j=(xi-1,xi)×(yj-1/2,jj+1/2) and then applying midpoint quadrature rule, we can obtain(22)-Δx+vi-1/2,jy-Δx-vi-1/2,jyxi-xi-1-Δy+vi-1/2,jy-Δy-vi-1/2,jyyj+1/2-yj-1/2+Δy-pi-1/2,j+1/2=fi-1/2,jy.By integrating (16) over cell Ωi-1/2,j+1/2=(xi-xi-1)×(yj,jj+1), we obtain its discrete form(23)Δx-vi,j+1/2x+Δy+vi-1/2,jy=0.Eqns. (21)–(23) form a linear equation system for the pressure and velocity, and the equation system defines the MAC scheme for the Stokes problem of (15)–(17). Note that the no-flow boundary condition supplies the zero velocity component that acts normal to any fluid-solid interface, and periodic conditions are applied along the domain’s external boundaries.We apply this numerical scheme to solve the set of rescaled Stokes equations in Eqns. (10)–(12). Our computer-generated media samples are 2-dimensional (2D); thus the set is solved with different forcings δij where i=1,2 and j=1,2. Upon obtaining the rescaled pore space velocity Nijy in the Y-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,(24)k=k11k12k21k22,inx⟹k∗=kmax⁡00kt,inx∗.


### 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 τh 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 x direction is(25)τh,x=vmagΩvxΩ,where 〈·〉Ω is the spatial average in the unit cell domain Ω, vmag is the magnitude of the fluid velocities (i.e., the fluid speed), and vx is the velocity component in the x direction. The spatial average of vx in Ω is(26)vxΩ=1V∫VfvxfdVf=1V∑i,jvxfi,jΔVf,where V is the volume of Ω (pore and solid space combined), Vf is the fluid space, and the integral is taken over the fluid space only. The spatial average for vmag is taken in the same way. The tortuosity in the y direction is computed in a similar manner as follows:(27)τh,y=vmagΩvyΩ,where the velocity is the solution for flow driven in the y direction. We note that (25) and (27) correspond to the hydraulic tortuosity as we have defined it in Section 2.2; that is τh=Le/Ls. If the term τh2 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 vx and vy explicitly. Instead, during the permeability computation for our 2D pore geometries, we obtain the rescaled pore space velocity fields Nij in the Y-cell, where i=1,2 and j=1,2. It can be seen that N11 and N21 are related to vx and vy, respectively, for flow that is driven by g=(g^,0) and that N12 and N22 are related to vx and vy, respectively, for flow that is driven by g=(0,g^). Thus, it is mathematically equivalent to compute the hydraulic tortuosities by(28)τh,x=N11+N21ΩN11Ω,τh,y=N12+N22ΩN22Ω,where Ω refers to the domain of the Y-cell, and the spatial averages of Nij are computed in the same way as shown by (26). Since the numerical scheme computes Nij on the fluid edge-centers, we compute the fluid cell-center quantities by averaging across the grid cells, where Nij=0 on the fluid-solid interfaces as per the boundary condition in (12).As mentioned, we diagonalize the permeability tensor to obtain the principal permeability components and the principal directionsθ and θ+π/2≔θ^. 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, τh in the principal direction θ is(29)τh,θ=NiθΩNθΩ,where Niθ=(N1θ)2+(N2θ)2 and where N1θ, N2θ, and Nθ are obtained by vector transformations from Nij. τh in the other principal direction θ^ is similarly obtained. The notation we use in our results section is τh,max⁡=τh,θ and τh,t=τh,θ^.


### 3.3. Effective Diffusion Coefficient and Diffusive Tortuosity

The problem is similarly defined by the two scales illustrated in Figure1. Within the pore space, the diffusive and convective transport of a species of concentration c is given by(30)∂∂tcε+1εvε,i∂∂xicε=∂∂xiDij∂∂xjcε,niDij∂∂xjcε=0,vε,i=0,where vε,i is divergence free (5), ε implies the variables are functions of the microscale, and Dij 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(31)cε=cx,xε=c0x,y+εc1x,y+ε2c2x,y+⋯and (8) are substituted into the governing transport equations (30). Upon collecting like order terms, the second-order term of c is(32)c1x,y=Pkx,y∂∂xkc0x, where c(0)(x) is the macroscale concentration and Pk(x,y) satisfies(33)∂∂yi∂∂yiPkx,y=0,(34)ni∂∂yiPkx,y=-nk,where (34) gives the boundary condition and ni 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, Dij=Dδij. Upon further algebra and averaging over the Y-cell, the macroscale transport equation is obtained as follows:(35)ϕ∂∂tcx=∂∂xiDikeff∗∂∂xkcx,where the effective diffusion coefficient Dikeff∗ is given by(36)Dikeff∗D=1Y∫Yδik+∂∂yiPkx,ydY.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(37)Dikeff∗=Dϕτd,ik,where τd,ik 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 by(38)τd,ik=Iik+1Yf∫Yf∂∂yiPkx,ydyf-1,where Yf=ϕY and where Pk(x,y) 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 Section3.1 to obtain the components τd,max⁡ and τd,t 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 τd∗, as will be shown by one of the following pore geometry examples (randomly distributed squares).


## 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 permeabilitykij is defined as the spatial average of arescaled pore space velocity Nij(y). 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 Figure1, 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:(4)-∂∂xipε+μ∂∂xj∂∂xjvε,i=-ρgj,(5)∂∂xivε,i=0,where i=1,2,3 and j=1,2,3 and repeated indices imply summation, p is pressure, μ is fluid viscosity, vi is the pore space fluid velocity, ρ is fluid density, and gj is gravitational acceleration. The boundary condition at the fluid-solid interface is given by(6)vi=0,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].Figure 1
Macroscale (a) and microscale (b) domains used in the homogenization problem. The solid structure in theY-cell illustrates the geometry of our first example: in-line array of solid shapes (i.e., circles or squares).In the Stokes equations,pε and vε,i can vary within the pore space and thus are functions of the macroscale (x) and microscale (ε); however μ, g, and ρ are treated as constants and are thus functions of the macroscale (x) only. By defining a fast variable y=x/ε, the cell of unit length ε is scaled into the Y-cell (shown in Figure 1(b)), where 0≤y≤1. A multiscale expansion is used to write pε and vε,i in terms of their leading and higher order terms as follows:(7)pε=px,xε=p0x,y+εp1x,y+⋯,(8)vε,i=vix,xε=vi0x,y+εvi1x,y+ε2vi2x,y+⋯.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(9)p1x,y=Qjy-ρgj+∂∂xjp0x,v2x,y=Nijyμ-ρgj+∂∂xjp0x,such that Qj(y) and Nij(y) satisfy(10)-∂∂yiQjy+∂∂yk∂∂ykNijy=δij,(11)∂∂yiNijy=0,with the fluid-solid boundary condition given by(12)Nijy=0.In the above formulation, Qj(y) and Nij(y) are functions of the microscale only and are independent of macroscale quantities μ, ρ, and gj. 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 Y-cell. In (10), δij is the Kronecker delta function (δij=1 for i=j, and δij=0 for i≠j).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 ofv(2)(x,y) in the cell is found to be divergence free; ∇·v2x,yy=0. This implies that the macroscopic flow equation is satisfied by Darcy’s equation (1), that is,(13)∂∂xikijμ∂∂xjpx-ρgj=0⟶∂∂xiuix=0,where p(x) is the macroscale pressure and ui(x) is the macroscale velocity, when the full permeability tensor for an n-dimensional problem is defined as(14)kij≔Nijyy=1Y∫YNijydY,where Y is the volume of the Y-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 Y-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=(vx,vy)T and the pressure p such that(15)-μΔv+∇p=f,inΩ≔0,Lx×0,Ly,(16)∇·v=0,inΩ,(17)v=0,on⁡∂Ω,where the velocity has a unique solution while the pressure is unique up to an additive constant. Here Lx and Ly denote the length of the rectangular domain Ω in the x- and y-directions, respectively.To construct a staggered grid finite difference scheme for the Stokes problem, we introduce a possibly nonuniform grid ofΩ as follows: (18)0=x0<x1<⋯<xM=Lx,0=y0<y1<⋯<yN=Ly.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 (xi+1/2,yj+1/2) is p(xi+1/2,yj+1/2) approximated by pi+1/2,j+1/2. Similarly, we approximate the x-velocity vx at each x-edge (fluid) center vx(xi,yj+1/2) by vi,j+1/2x, and the y-velocity vy at each y-edge (fluid) center vy(xi-1/2,yj) by vi-1/2,jy.By integrating thex-component of (15) over cell Ωi,j+1/2=(xi-1/2,xi+1/2)×(yj,yj+1), we have(19)-∫yjyj+1∂vx∂xxi+1/2,y-∂vx∂xxi-1/2,ydy+∫xi-1/2xi+1/2∂vx∂yx,yj+1-∂vx∂yx,yjdx+∫yjyj+1pxi+1/2,y-pxi-1/2,ydy=∫Ωi,j+1/2fxx,ydxdy.We then approximate those edge integrals by midpoint values and the volume integral by its center value and divide both sides byΩi,j+1/2=(xi+1/2-xi-1/2)(yj+1-yj).By introducing the following forward and backward difference operators:(20)Δx+vi,j+1/2x=vi+1,j+1/2x-vi,j+1/2xxi+1-xi,Δx-pi+1/2,j+1/2=pi+1/2,j+1/2-pi-1/2,j+1/2xi+1/2-xi-1/2and similar notations for operations in the y-direction, we can express the results as(21)-Δx+vi,j+1/2x-Δx-vi,j+1/2xxi+1/2-xi-1/2-Δy+vi,j+1/2x-Δy-vi,j+1/2xyj+1-yj+Δx-pi+1/2,j+1/2=fi,j+1/2x.Similarly, by integrating they-component of (15) over cell Ωi-1/2,j=(xi-1,xi)×(yj-1/2,jj+1/2) and then applying midpoint quadrature rule, we can obtain(22)-Δx+vi-1/2,jy-Δx-vi-1/2,jyxi-xi-1-Δy+vi-1/2,jy-Δy-vi-1/2,jyyj+1/2-yj-1/2+Δy-pi-1/2,j+1/2=fi-1/2,jy.By integrating (16) over cell Ωi-1/2,j+1/2=(xi-xi-1)×(yj,jj+1), we obtain its discrete form(23)Δx-vi,j+1/2x+Δy+vi-1/2,jy=0.Eqns. (21)–(23) form a linear equation system for the pressure and velocity, and the equation system defines the MAC scheme for the Stokes problem of (15)–(17). Note that the no-flow boundary condition supplies the zero velocity component that acts normal to any fluid-solid interface, and periodic conditions are applied along the domain’s external boundaries.We apply this numerical scheme to solve the set of rescaled Stokes equations in Eqns. (10)–(12). Our computer-generated media samples are 2-dimensional (2D); thus the set is solved with different forcings δij where i=1,2 and j=1,2. Upon obtaining the rescaled pore space velocity Nijy in the Y-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,(24)k=k11k12k21k22,inx⟹k∗=kmax⁡00kt,inx∗.


## 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 τh 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 x direction is(25)τh,x=vmagΩvxΩ,where 〈·〉Ω is the spatial average in the unit cell domain Ω, vmag is the magnitude of the fluid velocities (i.e., the fluid speed), and vx is the velocity component in the x direction. The spatial average of vx in Ω is(26)vxΩ=1V∫VfvxfdVf=1V∑i,jvxfi,jΔVf,where V is the volume of Ω (pore and solid space combined), Vf is the fluid space, and the integral is taken over the fluid space only. The spatial average for vmag is taken in the same way. The tortuosity in the y direction is computed in a similar manner as follows:(27)τh,y=vmagΩvyΩ,where the velocity is the solution for flow driven in the y direction. We note that (25) and (27) correspond to the hydraulic tortuosity as we have defined it in Section 2.2; that is τh=Le/Ls. If the term τh2 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 vx and vy explicitly. Instead, during the permeability computation for our 2D pore geometries, we obtain the rescaled pore space velocity fields Nij in the Y-cell, where i=1,2 and j=1,2. It can be seen that N11 and N21 are related to vx and vy, respectively, for flow that is driven by g=(g^,0) and that N12 and N22 are related to vx and vy, respectively, for flow that is driven by g=(0,g^). Thus, it is mathematically equivalent to compute the hydraulic tortuosities by(28)τh,x=N11+N21ΩN11Ω,τh,y=N12+N22ΩN22Ω,where Ω refers to the domain of the Y-cell, and the spatial averages of Nij are computed in the same way as shown by (26). Since the numerical scheme computes Nij on the fluid edge-centers, we compute the fluid cell-center quantities by averaging across the grid cells, where Nij=0 on the fluid-solid interfaces as per the boundary condition in (12).As mentioned, we diagonalize the permeability tensor to obtain the principal permeability components and the principal directionsθ and θ+π/2≔θ^. 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, τh in the principal direction θ is(29)τh,θ=NiθΩNθΩ,where Niθ=(N1θ)2+(N2θ)2 and where N1θ, N2θ, and Nθ are obtained by vector transformations from Nij. τh in the other principal direction θ^ is similarly obtained. The notation we use in our results section is τh,max⁡=τh,θ and τh,t=τh,θ^.


## 3.3. Effective Diffusion Coefficient and Diffusive Tortuosity

The problem is similarly defined by the two scales illustrated in Figure1. Within the pore space, the diffusive and convective transport of a species of concentration c is given by(30)∂∂tcε+1εvε,i∂∂xicε=∂∂xiDij∂∂xjcε,niDij∂∂xjcε=0,vε,i=0,where vε,i is divergence free (5), ε implies the variables are functions of the microscale, and Dij 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(31)cε=cx,xε=c0x,y+εc1x,y+ε2c2x,y+⋯and (8) are substituted into the governing transport equations (30). Upon collecting like order terms, the second-order term of c is(32)c1x,y=Pkx,y∂∂xkc0x, where c(0)(x) is the macroscale concentration and Pk(x,y) satisfies(33)∂∂yi∂∂yiPkx,y=0,(34)ni∂∂yiPkx,y=-nk,where (34) gives the boundary condition and ni 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, Dij=Dδij. Upon further algebra and averaging over the Y-cell, the macroscale transport equation is obtained as follows:(35)ϕ∂∂tcx=∂∂xiDikeff∗∂∂xkcx,where the effective diffusion coefficient Dikeff∗ is given by(36)Dikeff∗D=1Y∫Yδik+∂∂yiPkx,ydY.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(37)Dikeff∗=Dϕτd,ik,where τd,ik 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 by(38)τd,ik=Iik+1Yf∫Yf∂∂yiPkx,ydyf-1,where Yf=ϕY and where Pk(x,y) 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 Section3.1 to obtain the components τd,max⁡ and τd,t 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 τd∗, 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 (ϕ<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 Sections3.1, 3.2, and 3.3, respectively. Each pore geometry is generated by defining certain input parameters (i.e., radius of solid circle, side length of square, and density of randomly distributed squares), and these inputs are summarized in Table 4. By using a range of values for the input parameters, we obtain unit cells which we treat as representative elementary volumes over a range of porosities. After computing the permeability and tortuosity, we plot the data to observe permeability-porosity and tortuosity-porosity trends.Table 4
Input parameters used to generate different pore-geometries (or differentY-cells) and their corresponding porosity ranges. In this table, cell refers to a cell in the mesh used to discretize the Y-cell.
Geometry Initial mesh dimensions Inputs Porosity range In-line array of shapes (see Figure1) 50 × 50 cells (Y-cell length: 1 unit) (a) Circle diameter: 44 to 8 cells(b) Square length: 48 to 2 cells 0.33 < ϕ < 1 0.07<ϕ<1  Staggered array of squares (see Figure7) 100 × 50 cells (Y-cell length: Lb = 0.5 units, a/b = 1, La/Lb = 1 as per Figure 8) Square length: 40 to 8 cells 0.36 ≤ ϕ < 1  Randomly distributed squares (see Figure14) 100 × 100 cells (Y-cell length: 1 unit) Increasing density of squares (with constant square length: 0.01 unit = 1 cell) 0.45 ≤ ϕ < 1


### 4.1. In-Line Array of Uniform Shapes

The unit cell, orY-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 for0.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., kxx=kyy) 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 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.Figure 2
Permeability-porosity trend for in-line array of circles (i.e., infinitely long cylinders). Permeability units are nondimensional (i.e., unit length× unit length). The permeability is isotropic and is normalized by r2, where r is the solid circle radius. Gebart’s [47] analytical solution was given in (A.2).We obtained the hydraulic tortuosity in thex and y 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 x and y directions were computed to be τh,x=τh,y=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.Figure 3
Fluid flow past a solid circle,ϕ=0.71, driven by gx. 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 gy are same as shown but rotated 90 degrees.


#### 4.1.2. Diffusive Tortuosity

We obtained the full diffusive tortuosity tensor on increasingly refined meshes until the tortuosity converged within a0.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 τd 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 Pk (k=1,2 or x,y) for a unit cell of ϕ=0.71, which were obtained by solving the homogenized problem described in Section 3.3. The resulting diffusive tortuosity tensor is(39)τd=τd,xxτd,xyτd,yxτd,yy=1.2887-1×10-16-1×10-161.2887,where τd=τd′-1. Notice the tensor is symmetric and diagonal (i.e., τd,xy=τd,yx=0 to machine precision) and isotropic (i.e., τd,xx=τd,yy). 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 ∇2C=0 with Dirichlet boundary conditions on the right and left boundaries such that the concentration gradient across the unit cell is one.Figure 4
Number of grid cells required to reach convergedτd (using a convergence criteria of 0.1% relative error) for select porosities of in-line array of circles.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 Pk, and (c) magnitude of concentration gradient (∇xCx)2+(∇yCx)2 (diffusion speed) and direction. The fields with respect to y are the same as shown above but rotated by 90 degrees.
(a)
P x (b)
C x = P x + x s (c)
Diffusion speed and direction


#### 4.1.3. Comparison of Tortuosities for In-Line Array Unit Cells

The comparison between our isotropic results forτh and τd illustrated in Figure 6 allows us to emphasize that these tortuosity types are quantitatively different. This may be surprising at first, since the streamlines in Figures 3 and 5(c) appear to indicate that the flow fields are relatively analogous. However, the fluid speeds shown in Figure 3 are such that the fastest moving flow exists along the nontortuous pathway, while the fastest passive diffusive transport (or highest concentration gradient) shown in Figure 5(c) exists along the circle’s rounded fluid-solid boundary. This difference in speeds (or streamline weights) leads to a quantitative difference between these two tortuosity types.Figure 6
Tortuosity-porosity trends for in-line array of uniform shapes (circle or square). (a) Hydraulic tortuosity is a very weak function of porosity. (b) Diffusive tortuosity for in-line array of squares follows the analytical solution (τ=2-ϕ) from [48] closely, and both in-line array of squares and in-line array of circles match closely results by F. J. Valdés-Parada (provided via personal communication, December 3, 2016).
(a) (b)To further make our point, we quantitatively compare the hydraulic and diffusive streamlines in Figures3 and 5(c). The length λ(y) of 12 streamlines with starting positions given by (0,ys,j) is tabulated in Table 5. These lengths are computed by λj=lj/L, where lj is the length of the jth streamline and L=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.Table 5
Hydraulic and diffusive streamline comparison: in-line array of circles,ϕ=0.71. Flow is from left to right. (0,ys,j) indicate the starting location of the jth streamline.
Hydraulic streamlines Diffusive streamlines ( 0 , y s , j ) λ ( y j ) in Figure 3 λ ( y j ) in Figure 5(c)  0.05 1.0011 1.0010 0.15 1.0107 1.0103 0.25 1.0348 1.0348 0.35 1.0827 1.0852 0.45 1.1665 1.1814 0.495 1.2444 1.2732 0.505 1.2444 1.2732 0.55 1.1665 1.1814 0.65 1.0827 1.0852 0.75 1.0348 1.0348 0.85 1.0107 1.0103 0.95 1.0011 1.0010In Figure6, 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 (ϕc=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.


### 4.2. Staggered-Array of Uniform Shapes

Figure7 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, a/b=1 and La/Lb=1. By varying the solid square length, a range of porosities are obtained.Figure 7
Unit cell for staggered-array of uniform shapes (i.e., squares), represented by cells with periodic boundaries and porosityϕ. For all staggered-array examples: a/b=1 and La/Lb=1 (refer to Figure 8 for dimension details).
(a)
ϕ = 0.36 (b)
ϕ = 0.64 (c)
ϕ = 0.84Figure 8
A staggered (anisotropic) unit cell with dimensionsϕ=0.6, a/b=3, and La/Lb=1.


#### 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 of100×50, the meshes required refinement up to a maximum of 700×350 in order to meet the convergence criteria of krel  error<1%.This pore geometry configuration is anisotropic (i.e.,kxx≠kyy), 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:(40)γ≔kyykxx,where kyy>kxx for all porosities studied, as shown in Figure 10.Figure 9
Permeability-porosity trend for staggered-array of squares.Figure 10
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 thex direction than in the y direction. The hydraulic flow is more permeable in the y direction.The hydraulic tortuosity in thex and y 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 τh,x=1.3539 and τh,y=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 x direction. However, there is a straight flow path when the flow is driven in the y direction. Due to this characteristic, we observe that the hydraulic flow is more tortuous in the x direction than in the y direction, which is confirmed by the anisotropic ratio for hydraulic tortuosity plotted in Figure 10. By observing Figure 11, we notice that the fluid speeds are highest in the pore regions located between two square corners. We also notice in Figure 11(a) that the fluid speed is close to being constant along a sinusoidal flow path which appears to occupy most of the pore space, except towards the fluid-solid boundaries. The horizontally driven flow is interrupted by the center square and causes the flow path to diverge. However, in Figure 11(b), the vertically driven flow is not interrupted to the extent that would cause the flow path to become sinusoidal. Instead, the geometry of the pore space allows for a straight flow path to develop. We do not notice this sinusoidal flow path of constant speed, rather we notice that the speed is more localized (i.e., speed is highest between adjacent solid corners and lowest along solid square sides).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.
(a)
Flow inx direction (b)
Flow iny direction


#### 4.2.2. Diffusive Tortuosity

We obtained the full diffusive tortuosity tensor for this staggered geometry using a convergence criteria ofτd,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(41)τd=τd,xxτd,xyτd,yxτd,yy=1.3514-1×10-151×10-171.3092,where the tensor is symmetric and diagonal (i.e., τd,xy=τd,yx=0 to machine precision) and anisotropic (i.e., τd,xx≠τd,yy). In this geometry of staggered squares, τd,xx>τd,yy, and thus passive diffusive transport is more tortuous in the x direction due to an interrupted flow path.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 Pk (assuming concentration gradient across unit cell is one), and (c, f) magnitude of concentration gradient (diffusion speed) and direction.
(a)
P x (b)
C x = P x + x s (c)
Diffusion speed and direction (d)
P y (e)
C y = P y + y s (f)
Diffusion speed and direction


#### 4.2.3. Comparison of Tortuosities for Staggered-Array Unit Cells

Our results for hydraulic and diffusive tortuosity over the range of porosities are plotted in Figure13. 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 a/b=1, La/Lb=1), we do not have an explanation as to how τd could be smaller than 1, other than a possible typographical error. Since they reported values for the parameter ϕDeff/D rather than τd explicitly, it may not have been easily apparent to them that one of their values corresponded to τd<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 τd,xx 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 τd,yy 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 x 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, τh,x>τh,y for all porosities as shown in Figure 13, and τh,y is a very weak function of porosity. The anisotropic ratios of the tortuosities are(42)τh,yτh,x,τd,yyτd,xx 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 x direction than in the y direction. The hydraulic tortuosity shows a greater degree of anisotropy than diffusive tortuosity because τh,x≈1 for all porosities while τh,y increases towards the lower porosities.Table 6
Data comparison for diffusive tortuosity.
Unit cell specs. Kim et al.’s [43] results Our results ϕ a / b L a / L b ϕ D x x eff / D ϕ D y y eff / D ϕ τ d , x x ′ ϕ τ d , y y ′ 0.74 15 1 0.741 0.028 0.7396 0.0234 0.80 15 1 0.787 0.116 0.7959 0.1107 0.88 15 1 0.870 0.225 0.8675 0.2197 0.91 15 1 0.902 0.301 0.9112 0.3022  0.6 1 1 0.427 0.412 0.3863 0.4125 0.65 1 1 0.478 0.464 0.4397 0.4630 0.7 1 1 0.535 0.519 0.4940 0.5140 0.75 1 1 0.596 0.581 0.5625 0.5782 0.8 1 1 0.664 0.648 0.6439 0.6547 0.85 1 1 0.738 0.722 0.7221 0.7289 0.9 1 1 0.819 0.805 0.8064 0.8097 0.95 1 1 0.907 0.897 0.8976 0.8985  0.6 3 1 0.494 0.161 0.4814 0.1620 0.65 3 1 0.562 0.209 0.5465 0.2085 0.7 3 1 0.627 0.270 0.6187 0.2693 0.75 3 1 0.689 0.345 0.6747 0.3453 0.8 3 1 0.748 0.434 0.7350 0.4341 0.85 3 1 0.807 0.541 0.8012 0.5405 0.9 3 1 0.868 0.667 0.8627 0.6637 0.95 3 1 0.933 0.817 0.9282 0.8201Figure 13
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 they direction is a weak function of porosity.
(a)
Tortuosity in horizontal direction (b)
Tortuosity in vertical direction


### 4.3. Randomly Distributed Squares

Our third pore geometry generated comprised of randomly distributed squares, which is a type commonly used in literature (e.g., [13, 19, 26]). In those works, squares are freely overlapping; however, in our work, our square size is equivalent to the cell size of the initial mesh resolution, and thus no squares overlap. In this geometry, periodic boundary conditions were implemented by checking for the existence of any solid-fluid interfaces along the external boundaries of the domain, and, if detected, no-flow boundary conditions were applied at the interface (i.e., the velocity component acting normal to the fluid-solid interface was assigned a value of zero in the system of equations). We imposed the constraint given by Koza et al. [34] (i.e., the ratio of obstacle length to domain length < 0.01) which was recommended to avoid anisotropic-related statistical errors. This same constraint was also respected in Duda et al. [13]. Our geometries ranged from (effective) porosities of 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.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.
(a)
ϕ = 0 . 53 (b)
ϕ = 0 . 61 (c)
ϕ = 0 . 70 (d)
ϕ = 0 . 80 (e)
ϕ = 0 . 90


#### 4.3.1. Permeability and Hydraulic Tortuosity

We computed permeability using a convergence criteria ofkrel  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 systemx) 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).Figure 15
Permeability-porosity trend for randomly distributed squares. The full permeability tensor contains nonzero off-diagonal components (absolute values are plotted; however negative and positive quantities are indicated by different symbols) and is diagonalized to obtain the principal permeability components,kmax⁡ and kt. x and x∗ are related by rotation angle θk, which is plotted in Figure 16(a). The anisotropic ratio of the principal components is plotted in Figure 16(b).
(a)
Full tensor components, inx (b)
Diagonalized tensor components, inx∗Figure 16
Anisotropic ratios in geometries comprised of randomly distributed squares: (a) rotation to obtain diagonal tensors (k∗ and τd∗) is independent of porosity; (b) anisotropic ratios of diagonal tensors and anisotropic ratio of hydraulic tortuosity. At lower porosities (ϕ<0.6), anisotropic properties have been induced.
(a)
Rotation angle,θ (b)
Anisotropic ratiosWe 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 τh,x=1.3727 and τh,y=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 vmag>5×10-5) exist in the geometry, while the rest of the pore space is characterized by moderate (1×10-5<vmag<5×10-5) to slow (vmag<1×10-5) 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 y direction only, similar characteristics were found for flow in the x direction; thus we have omitted those illustrations. The computed hydraulic tortuosity at this porosity is essentially isotropic; that is, τh,x≈τh,y. 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.Figure 17
Flow (driven bygy) 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).
(a)
Fluid speed (b)
Close-up of boxed region in (a)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).


#### 4.3.2. Fit to Kozeny-Carman Equation

As presented in Section2.2, the Kozeny-Carman equation relates permeability k to porosity ϕ, hydraulic tortuosity τh, specific surface area S, and a fitting coefficient βk (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 βk 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 βk 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].Figure 19
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.Figure 20
Shape factor,βk, used to fit simulated permeability and hydraulic tortuosity data to Kozeny-Carman equation for geometries comprised of randomly distributed squares. βk appears to be constant between 0.7<ϕ<0.9.


#### 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τ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(43)τd=τd,xxτd,xyτd,yxτd,yy=2.7584-0.2914-0.29142.6175.Notice the tensor is symmetric (i.e., τd,xy=τd,yx) butnot 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.Figure 21
Fields for geometry of randomly distributed squares,ϕ=0.7: (a) solution to the homogenization problem, (b) concentration field transformed from Pk (assuming concentration gradient across unit cell is one), and (c, d) magnitude of concentration gradient (diffusion speed).
(a)
P x (b)
C x = P x + x s (c)
( ∇ x C x ) 2 + ( ∇ y C x ) 2 (d)
Close-up of boxed region in (c)Figure 22
Diffusive tortuosity-porosity trend for randomly distributed squares. The full diffusive tortuosity tensor contains nonzero off-diagonal components and is diagonalized to obtain the principal diffusive tortuosity components,τd,max⁡ and τd,t. x and {x∗} are related by rotation angle θd, which is plotted in Figure 16(a). The anisotropic ratio of the principal components is plotted in Figure 16(b).
(a)
Full tensor components, inx (b)
Diagonalized tensor components, in{x∗}


#### 4.3.4. Comparison of Tortuosities for Randomly Distributed Squares

Our hydraulic and diffusive tortuosity data are plotted in Figure23 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.Figure 23
Comparison of hydraulic (τh,x and τh,y) and diffusive (τd,xx and τd,yy) tortuosity trends for pore structures of randomly distributed squares, against trends reported in literature.
(a)
Tortuosity-porosity trends (b)
Close-up of data


#### 4.3.5. Induced Anisotropic Properties at Lower Porosity Geometries

Due to the method we used to generate the randomly distributed squares geometry, anisotropic properties were introduced into the pore structure, especially at lower porosities. This is evident in Figure24, 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 τd∗, 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).)Figure 24
Induced anisotropy in lower porosity geometries comprised of randomly distributed squares.
(a)
ϕ = 0.4525 (b)
ϕ = 0 . 7005Figure 25
Comparison of hydraulic and diffusive flows through a low porosity geometry comprised of randomly distributed squares:ϕ=0.45. (a, b) Fluid flow driven by gx. (c, d) Diffusive flow driven by ∇C=1 from west to east. “Hot spots” of fastest moving hydraulic flow are evident in channel centers, while diffusive speeds are fastest along fluid-solid boundaries.
(a)
Fluid speed (b)
Close-up of boxed region in (a) (c)
Diffusion speed (d)
Close-up of boxed region in (c)


## 4.1. In-Line Array of Uniform Shapes

The unit cell, orY-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 for0.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., kxx=kyy) 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 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.Figure 2
Permeability-porosity trend for in-line array of circles (i.e., infinitely long cylinders). Permeability units are nondimensional (i.e., unit length× unit length). The permeability is isotropic and is normalized by r2, where r is the solid circle radius. Gebart’s [47] analytical solution was given in (A.2).We obtained the hydraulic tortuosity in thex and y 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 x and y directions were computed to be τh,x=τh,y=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.Figure 3
Fluid flow past a solid circle,ϕ=0.71, driven by gx. 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 gy are same as shown but rotated 90 degrees.


### 4.1.2. Diffusive Tortuosity

We obtained the full diffusive tortuosity tensor on increasingly refined meshes until the tortuosity converged within a0.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 τd 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 Pk (k=1,2 or x,y) for a unit cell of ϕ=0.71, which were obtained by solving the homogenized problem described in Section 3.3. The resulting diffusive tortuosity tensor is(39)τd=τd,xxτd,xyτd,yxτd,yy=1.2887-1×10-16-1×10-161.2887,where τd=τd′-1. Notice the tensor is symmetric and diagonal (i.e., τd,xy=τd,yx=0 to machine precision) and isotropic (i.e., τd,xx=τd,yy). 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 ∇2C=0 with Dirichlet boundary conditions on the right and left boundaries such that the concentration gradient across the unit cell is one.Figure 4
Number of grid cells required to reach convergedτd (using a convergence criteria of 0.1% relative error) for select porosities of in-line array of circles.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 Pk, and (c) magnitude of concentration gradient (∇xCx)2+(∇yCx)2 (diffusion speed) and direction. The fields with respect to y are the same as shown above but rotated by 90 degrees.
(a)
P x (b)
C x = P x + x s (c)
Diffusion speed and direction


### 4.1.3. Comparison of Tortuosities for In-Line Array Unit Cells

The comparison between our isotropic results forτh and τd illustrated in Figure 6 allows us to emphasize that these tortuosity types are quantitatively different. This may be surprising at first, since the streamlines in Figures 3 and 5(c) appear to indicate that the flow fields are relatively analogous. However, the fluid speeds shown in Figure 3 are such that the fastest moving flow exists along the nontortuous pathway, while the fastest passive diffusive transport (or highest concentration gradient) shown in Figure 5(c) exists along the circle’s rounded fluid-solid boundary. This difference in speeds (or streamline weights) leads to a quantitative difference between these two tortuosity types.Figure 6
Tortuosity-porosity trends for in-line array of uniform shapes (circle or square). (a) Hydraulic tortuosity is a very weak function of porosity. (b) Diffusive tortuosity for in-line array of squares follows the analytical solution (τ=2-ϕ) from [48] closely, and both in-line array of squares and in-line array of circles match closely results by F. J. Valdés-Parada (provided via personal communication, December 3, 2016).
(a) (b)To further make our point, we quantitatively compare the hydraulic and diffusive streamlines in Figures3 and 5(c). The length λ(y) of 12 streamlines with starting positions given by (0,ys,j) is tabulated in Table 5. These lengths are computed by λj=lj/L, where lj is the length of the jth streamline and L=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.Table 5
Hydraulic and diffusive streamline comparison: in-line array of circles,ϕ=0.71. Flow is from left to right. (0,ys,j) indicate the starting location of the jth streamline.
Hydraulic streamlines Diffusive streamlines ( 0 , y s , j ) λ ( y j ) in Figure 3 λ ( y j ) in Figure 5(c)  0.05 1.0011 1.0010 0.15 1.0107 1.0103 0.25 1.0348 1.0348 0.35 1.0827 1.0852 0.45 1.1665 1.1814 0.495 1.2444 1.2732 0.505 1.2444 1.2732 0.55 1.1665 1.1814 0.65 1.0827 1.0852 0.75 1.0348 1.0348 0.85 1.0107 1.0103 0.95 1.0011 1.0010In Figure6, 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 (ϕc=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.


## 4.1.1. Permeability and Hydraulic Tortuosity

The full permeability tensor was computed for0.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., kxx=kyy) 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 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.Figure 2
Permeability-porosity trend for in-line array of circles (i.e., infinitely long cylinders). Permeability units are nondimensional (i.e., unit length× unit length). The permeability is isotropic and is normalized by r2, where r is the solid circle radius. Gebart’s [47] analytical solution was given in (A.2).We obtained the hydraulic tortuosity in thex and y 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 x and y directions were computed to be τh,x=τh,y=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.Figure 3
Fluid flow past a solid circle,ϕ=0.71, driven by gx. 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 gy are same as shown but rotated 90 degrees.


## 4.1.2. Diffusive Tortuosity

We obtained the full diffusive tortuosity tensor on increasingly refined meshes until the tortuosity converged within a0.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 τd 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 Pk (k=1,2 or x,y) for a unit cell of ϕ=0.71, which were obtained by solving the homogenized problem described in Section 3.3. The resulting diffusive tortuosity tensor is(39)τd=τd,xxτd,xyτd,yxτd,yy=1.2887-1×10-16-1×10-161.2887,where τd=τd′-1. Notice the tensor is symmetric and diagonal (i.e., τd,xy=τd,yx=0 to machine precision) and isotropic (i.e., τd,xx=τd,yy). 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 ∇2C=0 with Dirichlet boundary conditions on the right and left boundaries such that the concentration gradient across the unit cell is one.Figure 4
Number of grid cells required to reach convergedτd (using a convergence criteria of 0.1% relative error) for select porosities of in-line array of circles.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 Pk, and (c) magnitude of concentration gradient (∇xCx)2+(∇yCx)2 (diffusion speed) and direction. The fields with respect to y are the same as shown above but rotated by 90 degrees.
(a)
P x (b)
C x = P x + x s (c)
Diffusion speed and direction


## 4.1.3. Comparison of Tortuosities for In-Line Array Unit Cells

The comparison between our isotropic results forτh and τd illustrated in Figure 6 allows us to emphasize that these tortuosity types are quantitatively different. This may be surprising at first, since the streamlines in Figures 3 and 5(c) appear to indicate that the flow fields are relatively analogous. However, the fluid speeds shown in Figure 3 are such that the fastest moving flow exists along the nontortuous pathway, while the fastest passive diffusive transport (or highest concentration gradient) shown in Figure 5(c) exists along the circle’s rounded fluid-solid boundary. This difference in speeds (or streamline weights) leads to a quantitative difference between these two tortuosity types.Figure 6
Tortuosity-porosity trends for in-line array of uniform shapes (circle or square). (a) Hydraulic tortuosity is a very weak function of porosity. (b) Diffusive tortuosity for in-line array of squares follows the analytical solution (τ=2-ϕ) from [48] closely, and both in-line array of squares and in-line array of circles match closely results by F. J. Valdés-Parada (provided via personal communication, December 3, 2016).
(a) (b)To further make our point, we quantitatively compare the hydraulic and diffusive streamlines in Figures3 and 5(c). The length λ(y) of 12 streamlines with starting positions given by (0,ys,j) is tabulated in Table 5. These lengths are computed by λj=lj/L, where lj is the length of the jth streamline and L=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.Table 5
Hydraulic and diffusive streamline comparison: in-line array of circles,ϕ=0.71. Flow is from left to right. (0,ys,j) indicate the starting location of the jth streamline.
Hydraulic streamlines Diffusive streamlines ( 0 , y s , j ) λ ( y j ) in Figure 3 λ ( y j ) in Figure 5(c)  0.05 1.0011 1.0010 0.15 1.0107 1.0103 0.25 1.0348 1.0348 0.35 1.0827 1.0852 0.45 1.1665 1.1814 0.495 1.2444 1.2732 0.505 1.2444 1.2732 0.55 1.1665 1.1814 0.65 1.0827 1.0852 0.75 1.0348 1.0348 0.85 1.0107 1.0103 0.95 1.0011 1.0010In Figure6, 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 (ϕc=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.


## 4.2. Staggered-Array of Uniform Shapes

Figure7 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, a/b=1 and La/Lb=1. By varying the solid square length, a range of porosities are obtained.Figure 7
Unit cell for staggered-array of uniform shapes (i.e., squares), represented by cells with periodic boundaries and porosityϕ. For all staggered-array examples: a/b=1 and La/Lb=1 (refer to Figure 8 for dimension details).
(a)
ϕ = 0.36 (b)
ϕ = 0.64 (c)
ϕ = 0.84Figure 8
A staggered (anisotropic) unit cell with dimensionsϕ=0.6, a/b=3, and La/Lb=1.


### 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 of100×50, the meshes required refinement up to a maximum of 700×350 in order to meet the convergence criteria of krel  error<1%.This pore geometry configuration is anisotropic (i.e.,kxx≠kyy), 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:(40)γ≔kyykxx,where kyy>kxx for all porosities studied, as shown in Figure 10.Figure 9
Permeability-porosity trend for staggered-array of squares.Figure 10
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 thex direction than in the y direction. The hydraulic flow is more permeable in the y direction.The hydraulic tortuosity in thex and y 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 τh,x=1.3539 and τh,y=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 x direction. However, there is a straight flow path when the flow is driven in the y direction. Due to this characteristic, we observe that the hydraulic flow is more tortuous in the x direction than in the y direction, which is confirmed by the anisotropic ratio for hydraulic tortuosity plotted in Figure 10. By observing Figure 11, we notice that the fluid speeds are highest in the pore regions located between two square corners. We also notice in Figure 11(a) that the fluid speed is close to being constant along a sinusoidal flow path which appears to occupy most of the pore space, except towards the fluid-solid boundaries. The horizontally driven flow is interrupted by the center square and causes the flow path to diverge. However, in Figure 11(b), the vertically driven flow is not interrupted to the extent that would cause the flow path to become sinusoidal. Instead, the geometry of the pore space allows for a straight flow path to develop. We do not notice this sinusoidal flow path of constant speed, rather we notice that the speed is more localized (i.e., speed is highest between adjacent solid corners and lowest along solid square sides).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.
(a)
Flow inx direction (b)
Flow iny direction


### 4.2.2. Diffusive Tortuosity

We obtained the full diffusive tortuosity tensor for this staggered geometry using a convergence criteria ofτd,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(41)τd=τd,xxτd,xyτd,yxτd,yy=1.3514-1×10-151×10-171.3092,where the tensor is symmetric and diagonal (i.e., τd,xy=τd,yx=0 to machine precision) and anisotropic (i.e., τd,xx≠τd,yy). In this geometry of staggered squares, τd,xx>τd,yy, and thus passive diffusive transport is more tortuous in the x direction due to an interrupted flow path.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 Pk (assuming concentration gradient across unit cell is one), and (c, f) magnitude of concentration gradient (diffusion speed) and direction.
(a)
P x (b)
C x = P x + x s (c)
Diffusion speed and direction (d)
P y (e)
C y = P y + y s (f)
Diffusion speed and direction


### 4.2.3. Comparison of Tortuosities for Staggered-Array Unit Cells

Our results for hydraulic and diffusive tortuosity over the range of porosities are plotted in Figure13. 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 a/b=1, La/Lb=1), we do not have an explanation as to how τd could be smaller than 1, other than a possible typographical error. Since they reported values for the parameter ϕDeff/D rather than τd explicitly, it may not have been easily apparent to them that one of their values corresponded to τd<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 τd,xx 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 τd,yy 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 x 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, τh,x>τh,y for all porosities as shown in Figure 13, and τh,y is a very weak function of porosity. The anisotropic ratios of the tortuosities are(42)τh,yτh,x,τd,yyτd,xx 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 x direction than in the y direction. The hydraulic tortuosity shows a greater degree of anisotropy than diffusive tortuosity because τh,x≈1 for all porosities while τh,y increases towards the lower porosities.Table 6
Data comparison for diffusive tortuosity.
Unit cell specs. Kim et al.’s [43] results Our results ϕ a / b L a / L b ϕ D x x eff / D ϕ D y y eff / D ϕ τ d , x x ′ ϕ τ d , y y ′ 0.74 15 1 0.741 0.028 0.7396 0.0234 0.80 15 1 0.787 0.116 0.7959 0.1107 0.88 15 1 0.870 0.225 0.8675 0.2197 0.91 15 1 0.902 0.301 0.9112 0.3022  0.6 1 1 0.427 0.412 0.3863 0.4125 0.65 1 1 0.478 0.464 0.4397 0.4630 0.7 1 1 0.535 0.519 0.4940 0.5140 0.75 1 1 0.596 0.581 0.5625 0.5782 0.8 1 1 0.664 0.648 0.6439 0.6547 0.85 1 1 0.738 0.722 0.7221 0.7289 0.9 1 1 0.819 0.805 0.8064 0.8097 0.95 1 1 0.907 0.897 0.8976 0.8985  0.6 3 1 0.494 0.161 0.4814 0.1620 0.65 3 1 0.562 0.209 0.5465 0.2085 0.7 3 1 0.627 0.270 0.6187 0.2693 0.75 3 1 0.689 0.345 0.6747 0.3453 0.8 3 1 0.748 0.434 0.7350 0.4341 0.85 3 1 0.807 0.541 0.8012 0.5405 0.9 3 1 0.868 0.667 0.8627 0.6637 0.95 3 1 0.933 0.817 0.9282 0.8201Figure 13
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 they direction is a weak function of porosity.
(a)
Tortuosity in horizontal direction (b)
Tortuosity in vertical direction


## 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 of100×50, the meshes required refinement up to a maximum of 700×350 in order to meet the convergence criteria of krel  error<1%.This pore geometry configuration is anisotropic (i.e.,kxx≠kyy), 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:(40)γ≔kyykxx,where kyy>kxx for all porosities studied, as shown in Figure 10.Figure 9
Permeability-porosity trend for staggered-array of squares.Figure 10
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 thex direction than in the y direction. The hydraulic flow is more permeable in the y direction.The hydraulic tortuosity in thex and y 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 τh,x=1.3539 and τh,y=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 x direction. However, there is a straight flow path when the flow is driven in the y direction. Due to this characteristic, we observe that the hydraulic flow is more tortuous in the x direction than in the y direction, which is confirmed by the anisotropic ratio for hydraulic tortuosity plotted in Figure 10. By observing Figure 11, we notice that the fluid speeds are highest in the pore regions located between two square corners. We also notice in Figure 11(a) that the fluid speed is close to being constant along a sinusoidal flow path which appears to occupy most of the pore space, except towards the fluid-solid boundaries. The horizontally driven flow is interrupted by the center square and causes the flow path to diverge. However, in Figure 11(b), the vertically driven flow is not interrupted to the extent that would cause the flow path to become sinusoidal. Instead, the geometry of the pore space allows for a straight flow path to develop. We do not notice this sinusoidal flow path of constant speed, rather we notice that the speed is more localized (i.e., speed is highest between adjacent solid corners and lowest along solid square sides).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.
(a)
Flow inx direction (b)
Flow iny direction


## 4.2.2. Diffusive Tortuosity

We obtained the full diffusive tortuosity tensor for this staggered geometry using a convergence criteria ofτd,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(41)τd=τd,xxτd,xyτd,yxτd,yy=1.3514-1×10-151×10-171.3092,where the tensor is symmetric and diagonal (i.e., τd,xy=τd,yx=0 to machine precision) and anisotropic (i.e., τd,xx≠τd,yy). In this geometry of staggered squares, τd,xx>τd,yy, and thus passive diffusive transport is more tortuous in the x direction due to an interrupted flow path.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 Pk (assuming concentration gradient across unit cell is one), and (c, f) magnitude of concentration gradient (diffusion speed) and direction.
(a)
P x (b)
C x = P x + x s (c)
Diffusion speed and direction (d)
P y (e)
C y = P y + y s (f)
Diffusion speed and direction


## 4.2.3. Comparison of Tortuosities for Staggered-Array Unit Cells

Our results for hydraulic and diffusive tortuosity over the range of porosities are plotted in Figure13. 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 a/b=1, La/Lb=1), we do not have an explanation as to how τd could be smaller than 1, other than a possible typographical error. Since they reported values for the parameter ϕDeff/D rather than τd explicitly, it may not have been easily apparent to them that one of their values corresponded to τd<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 τd,xx 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 τd,yy 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 x 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, τh,x>τh,y for all porosities as shown in Figure 13, and τh,y is a very weak function of porosity. The anisotropic ratios of the tortuosities are(42)τh,yτh,x,τd,yyτd,xx 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 x direction than in the y direction. The hydraulic tortuosity shows a greater degree of anisotropy than diffusive tortuosity because τh,x≈1 for all porosities while τh,y increases towards the lower porosities.Table 6
Data comparison for diffusive tortuosity.
Unit cell specs. Kim et al.’s [43] results Our results ϕ a / b L a / L b ϕ D x x eff / D ϕ D y y eff / D ϕ τ d , x x ′ ϕ τ d , y y ′ 0.74 15 1 0.741 0.028 0.7396 0.0234 0.80 15 1 0.787 0.116 0.7959 0.1107 0.88 15 1 0.870 0.225 0.8675 0.2197 0.91 15 1 0.902 0.301 0.9112 0.3022  0.6 1 1 0.427 0.412 0.3863 0.4125 0.65 1 1 0.478 0.464 0.4397 0.4630 0.7 1 1 0.535 0.519 0.4940 0.5140 0.75 1 1 0.596 0.581 0.5625 0.5782 0.8 1 1 0.664 0.648 0.6439 0.6547 0.85 1 1 0.738 0.722 0.7221 0.7289 0.9 1 1 0.819 0.805 0.8064 0.8097 0.95 1 1 0.907 0.897 0.8976 0.8985  0.6 3 1 0.494 0.161 0.4814 0.1620 0.65 3 1 0.562 0.209 0.5465 0.2085 0.7 3 1 0.627 0.270 0.6187 0.2693 0.75 3 1 0.689 0.345 0.6747 0.3453 0.8 3 1 0.748 0.434 0.7350 0.4341 0.85 3 1 0.807 0.541 0.8012 0.5405 0.9 3 1 0.868 0.667 0.8627 0.6637 0.95 3 1 0.933 0.817 0.9282 0.8201Figure 13
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 they direction is a weak function of porosity.
(a)
Tortuosity in horizontal direction (b)
Tortuosity in vertical direction


## 4.3. Randomly Distributed Squares

Our third pore geometry generated comprised of randomly distributed squares, which is a type commonly used in literature (e.g., [13, 19, 26]). In those works, squares are freely overlapping; however, in our work, our square size is equivalent to the cell size of the initial mesh resolution, and thus no squares overlap. In this geometry, periodic boundary conditions were implemented by checking for the existence of any solid-fluid interfaces along the external boundaries of the domain, and, if detected, no-flow boundary conditions were applied at the interface (i.e., the velocity component acting normal to the fluid-solid interface was assigned a value of zero in the system of equations). We imposed the constraint given by Koza et al. [34] (i.e., the ratio of obstacle length to domain length < 0.01) which was recommended to avoid anisotropic-related statistical errors. This same constraint was also respected in Duda et al. [13]. Our geometries ranged from (effective) porosities of 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.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.
(a)
ϕ = 0 . 53 (b)
ϕ = 0 . 61 (c)
ϕ = 0 . 70 (d)
ϕ = 0 . 80 (e)
ϕ = 0 . 90


### 4.3.1. Permeability and Hydraulic Tortuosity

We computed permeability using a convergence criteria ofkrel  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 systemx) 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).Figure 15
Permeability-porosity trend for randomly distributed squares. The full permeability tensor contains nonzero off-diagonal components (absolute values are plotted; however negative and positive quantities are indicated by different symbols) and is diagonalized to obtain the principal permeability components,kmax⁡ and kt. x and x∗ are related by rotation angle θk, which is plotted in Figure 16(a). The anisotropic ratio of the principal components is plotted in Figure 16(b).
(a)
Full tensor components, inx (b)
Diagonalized tensor components, inx∗Figure 16
Anisotropic ratios in geometries comprised of randomly distributed squares: (a) rotation to obtain diagonal tensors (k∗ and τd∗) is independent of porosity; (b) anisotropic ratios of diagonal tensors and anisotropic ratio of hydraulic tortuosity. At lower porosities (ϕ<0.6), anisotropic properties have been induced.
(a)
Rotation angle,θ (b)
Anisotropic ratiosWe 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 τh,x=1.3727 and τh,y=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 vmag>5×10-5) exist in the geometry, while the rest of the pore space is characterized by moderate (1×10-5<vmag<5×10-5) to slow (vmag<1×10-5) 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 y direction only, similar characteristics were found for flow in the x direction; thus we have omitted those illustrations. The computed hydraulic tortuosity at this porosity is essentially isotropic; that is, τh,x≈τh,y. 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.Figure 17
Flow (driven bygy) 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).
(a)
Fluid speed (b)
Close-up of boxed region in (a)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).


### 4.3.2. Fit to Kozeny-Carman Equation

As presented in Section2.2, the Kozeny-Carman equation relates permeability k to porosity ϕ, hydraulic tortuosity τh, specific surface area S, and a fitting coefficient βk (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 βk 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 βk 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].Figure 19
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.Figure 20
Shape factor,βk, used to fit simulated permeability and hydraulic tortuosity data to Kozeny-Carman equation for geometries comprised of randomly distributed squares. βk appears to be constant between 0.7<ϕ<0.9.


### 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τ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(43)τd=τd,xxτd,xyτd,yxτd,yy=2.7584-0.2914-0.29142.6175.Notice the tensor is symmetric (i.e., τd,xy=τd,yx) butnot 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.Figure 21
Fields for geometry of randomly distributed squares,ϕ=0.7: (a) solution to the homogenization problem, (b) concentration field transformed from Pk (assuming concentration gradient across unit cell is one), and (c, d) magnitude of concentration gradient (diffusion speed).
(a)
P x (b)
C x = P x + x s (c)
( ∇ x C x ) 2 + ( ∇ y C x ) 2 (d)
Close-up of boxed region in (c)Figure 22
Diffusive tortuosity-porosity trend for randomly distributed squares. The full diffusive tortuosity tensor contains nonzero off-diagonal components and is diagonalized to obtain the principal diffusive tortuosity components,τd,max⁡ and τd,t. x and {x∗} are related by rotation angle θd, which is plotted in Figure 16(a). The anisotropic ratio of the principal components is plotted in Figure 16(b).
(a)
Full tensor components, inx (b)
Diagonalized tensor components, in{x∗}


### 4.3.4. Comparison of Tortuosities for Randomly Distributed Squares

Our hydraulic and diffusive tortuosity data are plotted in Figure23 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.Figure 23
Comparison of hydraulic (τh,x and τh,y) and diffusive (τd,xx and τd,yy) tortuosity trends for pore structures of randomly distributed squares, against trends reported in literature.
(a)
Tortuosity-porosity trends (b)
Close-up of data


### 4.3.5. Induced Anisotropic Properties at Lower Porosity Geometries

Due to the method we used to generate the randomly distributed squares geometry, anisotropic properties were introduced into the pore structure, especially at lower porosities. This is evident in Figure24, 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 τd∗, 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).)Figure 24
Induced anisotropy in lower porosity geometries comprised of randomly distributed squares.
(a)
ϕ = 0.4525 (b)
ϕ = 0 . 7005Figure 25
Comparison of hydraulic and diffusive flows through a low porosity geometry comprised of randomly distributed squares:ϕ=0.45. (a, b) Fluid flow driven by gx. (c, d) Diffusive flow driven by ∇C=1 from west to east. “Hot spots” of fastest moving hydraulic flow are evident in channel centers, while diffusive speeds are fastest along fluid-solid boundaries.
(a)
Fluid speed (b)
Close-up of boxed region in (a) (c)
Diffusion speed (d)
Close-up of boxed region in (c)


## 4.3.1. Permeability and Hydraulic Tortuosity

We computed permeability using a convergence criteria ofkrel  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 systemx) 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).Figure 15
Permeability-porosity trend for randomly distributed squares. The full permeability tensor contains nonzero off-diagonal components (absolute values are plotted; however negative and positive quantities are indicated by different symbols) and is diagonalized to obtain the principal permeability components,kmax⁡ and kt. x and x∗ are related by rotation angle θk, which is plotted in Figure 16(a). The anisotropic ratio of the principal components is plotted in Figure 16(b).
(a)
Full tensor components, inx (b)
Diagonalized tensor components, inx∗Figure 16
Anisotropic ratios in geometries comprised of randomly distributed squares: (a) rotation to obtain diagonal tensors (k∗ and τd∗) is independent of porosity; (b) anisotropic ratios of diagonal tensors and anisotropic ratio of hydraulic tortuosity. At lower porosities (ϕ<0.6), anisotropic properties have been induced.
(a)
Rotation angle,θ (b)
Anisotropic ratiosWe 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 τh,x=1.3727 and τh,y=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 vmag>5×10-5) exist in the geometry, while the rest of the pore space is characterized by moderate (1×10-5<vmag<5×10-5) to slow (vmag<1×10-5) 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 y direction only, similar characteristics were found for flow in the x direction; thus we have omitted those illustrations. The computed hydraulic tortuosity at this porosity is essentially isotropic; that is, τh,x≈τh,y. 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.Figure 17
Flow (driven bygy) 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).
(a)
Fluid speed (b)
Close-up of boxed region in (a)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).


## 4.3.2. Fit to Kozeny-Carman Equation

As presented in Section2.2, the Kozeny-Carman equation relates permeability k to porosity ϕ, hydraulic tortuosity τh, specific surface area S, and a fitting coefficient βk (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 βk 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 βk 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].Figure 19
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.Figure 20
Shape factor,βk, used to fit simulated permeability and hydraulic tortuosity data to Kozeny-Carman equation for geometries comprised of randomly distributed squares. βk appears to be constant between 0.7<ϕ<0.9.


## 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τ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(43)τd=τd,xxτd,xyτd,yxτd,yy=2.7584-0.2914-0.29142.6175.Notice the tensor is symmetric (i.e., τd,xy=τd,yx) butnot 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.Figure 21
Fields for geometry of randomly distributed squares,ϕ=0.7: (a) solution to the homogenization problem, (b) concentration field transformed from Pk (assuming concentration gradient across unit cell is one), and (c, d) magnitude of concentration gradient (diffusion speed).
(a)
P x (b)
C x = P x + x s (c)
( ∇ x C x ) 2 + ( ∇ y C x ) 2 (d)
Close-up of boxed region in (c)Figure 22
Diffusive tortuosity-porosity trend for randomly distributed squares. The full diffusive tortuosity tensor contains nonzero off-diagonal components and is diagonalized to obtain the principal diffusive tortuosity components,τd,max⁡ and τd,t. x and {x∗} are related by rotation angle θd, which is plotted in Figure 16(a). The anisotropic ratio of the principal components is plotted in Figure 16(b).
(a)
Full tensor components, inx (b)
Diagonalized tensor components, in{x∗}


## 4.3.4. Comparison of Tortuosities for Randomly Distributed Squares

Our hydraulic and diffusive tortuosity data are plotted in Figure23 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.Figure 23
Comparison of hydraulic (τh,x and τh,y) and diffusive (τd,xx and τd,yy) tortuosity trends for pore structures of randomly distributed squares, against trends reported in literature.
(a)
Tortuosity-porosity trends (b)
Close-up of data


## 4.3.5. Induced Anisotropic Properties at Lower Porosity Geometries

Due to the method we used to generate the randomly distributed squares geometry, anisotropic properties were introduced into the pore structure, especially at lower porosities. This is evident in Figure24, 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 τd∗, 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).)Figure 24
Induced anisotropy in lower porosity geometries comprised of randomly distributed squares.
(a)
ϕ = 0.4525 (b)
ϕ = 0 . 7005Figure 25
Comparison of hydraulic and diffusive flows through a low porosity geometry comprised of randomly distributed squares:ϕ=0.45. (a, b) Fluid flow driven by gx. (c, d) Diffusive flow driven by ∇C=1 from west to east. “Hot spots” of fastest moving hydraulic flow are evident in channel centers, while diffusive speeds are fastest along fluid-solid boundaries.
(a)
Fluid speed (b)
Close-up of boxed region in (a) (c)
Diffusion speed (d)
Close-up of boxed region in (c)


## 5. Conclusion

In this study, we computed the effective properties known as permeability, hydraulic tortuosity, and diffusive tortuosity, using well-known and commonly employed formulations based on the method of homogenization. While we noted that past work has recognized the difference between tortuosity types, we emphasized that a few studies have failed to distinguish the difference and have mistakenly used them interchangeably. In this work, hydraulic tortuosity was computed based on the approach given in Koponen et al. [19, 26]. Their approach uses pore space velocity fields v; however we used the rescaled velocity fields Nij which we obtained by solving a set of rescaled Stokes equations formulated by homogenization ((10)–(12)). Using Nij in (28) is mathematically equivalent to using v 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 βk.The main findings from the geometries studied are as follows:(1)
Hydraulic tortuosity isnot 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τd>τh 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 τd≤τh. 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 Figure10). 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 τd became anisotropic by varying degrees. It is interesting to note that τh 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 Figures10 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.



---

*Source: 4517259-2017-02-13.xml*