The Path Tracking Method as an alternative for tortuosity determination in granular beds

Tortuosity is one of the most elusive parameters of porous media. The fundamental question is whether it may be computed from the geometry only or the transport equations must be solved first. Elimination of the transport equations would significantly decrease the computation time and memory consumption and thus allow investigating larger samples. We compare the geometric to hydraulic tortuosity of a sphere-packed porous media. We applied the Discrete Element Method to generate a set of virtual beds based on experimental data taking into account the real porosity and particle distribution, the Lattice Boltzmann Method to compute the hydraulic tortuosity and geometrical approach, i.e. so-called Path Tracking Method, to calculate the geometrical tortuosity. Our study shows that the calculation time can be reduced from hours (if the LBM is used) to seconds (if the PTM is applied) without losing the accuracy of the final results. The relative error between average values of the tortuosity obtained for both used methods is less than 3%. We show that the applied geometrical method may serve as an attractive alternative to hydraulic tortuosity, particularly in large granular systems.


Introduction
The fluid flow through porous media is a complex, multiscale phenomena. Fluid particles penetrate pore space and draw paths of various shapes and lengths. If those paths are short, the drag of porous matrix is relatively low. If particles turn and travel on long and wiggly trajectories, the drag is high, what results in low permeability of the media. In order to take this effect into account, Kozeny [24] introduced a new physical quantity named tortuosity. He applied this parameter to correct the value of the hydraulic drop (which in turn depends on the path length) occurring during fluid flow through a porous body (Fig. 1a). The same parameter was then used by Carman [9] who used the tortuosity as a correction factor of the flow velocity in a porous medium  1b). Additionally, Carman proposed a new formula for calculation of the so-called specific surface. The Kozeny formula modified by the Carman is now widely known as the Kozeny-Carman equation, which is destined to predict the pressure drop in porous media: where p-pressure [Pa], x-coordinate along which the pressure drop occurs [m], C K C -Kozeny-Carman pore shape factor (a model constant), which should be equal to 5.0 (no unit) [9], T f -tortuosity factor (no unit) defined as the square of the tortuosity T (no unit), φ-porosity (no unit), μ-dynamic viscosity coefficient [kg/(m·s)], v f -filtration velocity [m/s]. Tortuosity T is defined as the ratio of the actual path length inside channels (L p ) to the thickness of the porous medium (L 0 ): In a general case, tortuosity may be understood as a geometrical quantity (geometric tortuosity, T g ) or as a flow Fig. 1 Visualisation of ideas proposed by: a Kozeny-the real path length inside pore channels (L p ) is higher than the thickness of the porous medium (L 0 ); b Carman-the real fluid velocity inside pore channels (v p ) is higher that the space averaged filtration velocity (v f ) (a) (b) property (hydraulic tortuosity, T h ). Other kinds of tortuosity are also known, e.g. diffusional or electric tortuosity. Important is that obtaining the tortuosity value in a specific case is still very difficult and for this reason the application of the Kozeny-Carman equation (as well as many other formulas where tortuosity occurs) is limited. Due to such a state, all studies related to new ways of obtaining this physical quantity are needed. It should be mentioned that in the literature analytical formulas used to calculate the tortuosity in granular beds based on porosity may be found. However, the range of results obtained from these formulas is quite large [44], what leads to difficulties in chosing appropriate formula in a specific case.
The hydraulic tortuosity may be calculated from a velocity field as follows: where v x [m/s] is the velocity component in direction of macroscopic flow in the porous material and v [m/s] is the velocity magnitude. The sum goes over pore space [13,22]. The velocity field needed for the application of the formula (3) may be obtained experimentally (e.g. by the use of the Particle Image Velocimetry) or numerically using Computational Fluid Dynamics. In the last case, Finite Difference Method, Finite Volume Method, Immersed Boundary Method, Lattice Boltzmann Method or even Finite Element Method may be applied. In the context of porous media, Lattice Boltzmann Method is particularly convenient due it the feature of easy mapping of complex geometries [29][30][31]38].
The methodology of calculating the hydraulic tortuosity on the basis of a velocity field was used before. It used directly velocity field [22] and path lines [32]. There, the hydraulic tortuosity was measured in 2D random generated pore structures consisting of squares and rectangles. Similar investigations were then performed by Nabovati and Sousa [35] and Duda et al. [13]. The results obtained by these Authors are quite consistent, but cannot be directly compared with results of this article due to differences in the porous media model. The same methodology was used by Aaltosalmi [3]. His work was focused on Lattice Boltzmann Method as such and on its application as for the simulation of fluid flows through materials such as paper, fabric, sandstone and others. Lane [26] performed a study on the hydraulic and the electric tortuosity in 2D sphere packs with a normal distribution.
Papers directly related to the investigations of the hydraulic tortuosity in granular porous media are rare; e.g. Wang [49] investigated the relation between the porosity and the tortuosity in a channel filled with spherical particles with a constant diameter. The number of particles was not stated, but it can be estimated at a few thousands. In these simulations (performed by the use of the Discrete Element Method) the particles were randomly added over the bed in a calculation loop and next they fall down by gravity. In the second stage, the BGK D3Q19 Lattice Boltzmann Model was applied. On the basis of figures shown there it may by stated that for the porosity equal to 0.4 (which value corresponds approximately to the porosity of a bed sample used in this article), the hydraulic tortuosity is approx. 1.1975. Data on computational time needed for performing simulations is not given.
Summarizing, papers related to the investigations of the hydraulic tortuosity in large porous beds created with the use of the Discrete Element Method, taking into account the particle size distribution, are not available. In the article, we propose such a methodology, however the relation between the particle size distribution and the tortuosity is not studied.
The geometric tortuosity may be obtained in three ways: (1) through correlations with other parameters characterizing the geometry of a porous body; (2) experimentally and (3) by using algorithms that track the direction or shape of the porous space. In the first approach it is usually assumed that the geometric tortuosity is a direct function of the porosity. It is the most popular approach to calculate the tortuosity, and many functions (empirical or derived analytically), in which the relationships between these quantities are proposed, may be found in literature [2,4,21,44,48]. In the experimental investigations, the computed tomography and image analysis [15,18,50], acoustic methods [19,20,27], optical methods [36] as well as other methods [17] are used. In case of the last mentioned possibility, firstly the geometry of the pore part must be defined (or determined), and then its arrangement in the space is tracked.
Yu and Li [52] analyzed a two-dimensional system of connected squares and proposed an analytical formula to calculate the tortuosity. Starley et al. [46] propose a three-dimensional tracer metric numerical model, based on mutually perpendicular cylinders, to estimate the tortuosity. Montes et al. [34] used the so-called Pathfinding A* algorithm to analyse electrical and thermal tortuosity in powder compacts (however, details concerning this algorithm were not given). In the literature, several other algoritms refering to the analysis the geometric tortuosity (or in fact: any path in a partially limited space) may be also found, however they are related to other research areas, such as medicine, geodesy, games theory, logistics, environmental science and other, and co they cannot be applied or even compared with the investigations described in the article. The numerical algorithm used here to calculate the geometric tortuosity in granular beds is called the Path Tracking Method (see Sect. 2.5) and was developed by Sobieski in 2009 [42,45] (details in Sect. 2.5).
In the paper, two different methods of determining the tortuosity with the use of the same virtual beds are comparedthese are the methodology proposed by Koponen et al. [22,23] and a methodology based on numerical analysis of the shapes of pore channels. Our investigations contain following steps: measurement of particle diameters and porosity in a sample of a real granular bed; specifying the distribution of diameters in this bed; division of the bed into fractions according to the specified distribution; generating the virtual bed representing the real bed sample (by simulations with the use of Discrete Element Method); obtaining the hydraulic tortuosity (by simulations with the use of Lattice Boltzmann Method); obtaining the geometric tortuosity (by the use of the above mentioned Path Tracking Method) and the comparison of results. All these steps are described below. We discuss not only the issues related directly to the tortuosity, but also with the porosity (as indicates the literature review it is the main factor impacting the tortuosity value) and with chosen features of the used methodologies. We show that the Path Tracking Method may be very useful due to the low demand for computing power (in opposition to Lattice Boltzmann Methods), relatively simple implementation (appropriate references are given) and the possibility of use in large granular beds. The last feature is particularly important, because the geometric method may be used even in cases, in which the application of the LBM approach is impossible.

Materials
In order to generate numerical samples of granular beds we need to know both particle sizes and the porosity of the sample. Both will be used as an input data for packing procedures. We used a real set of the SiLibeads Glass beads (Type S). First, we measured their diameters along two axes using the micrometer screw with ± 0.01 [mm] accuracy. We found that diameters of 100 beads follow Gaussian distribution at a significance level 0.01 of the Kolmogorov-Smirnov test [28]. The mean diameter of 100 beads was foundd = 6.072 [mm] with standard deviation σ = 0.051 [mm].
To test porosity of our sphere-packed media we loosely packed those glass beads into a graduated 200 [ml] cylinder and filled it slowly with water. We measured porosity by dividing the water volume (V 0 ) by the total volume of cylinder (V c ): φ = V 0 /V c . We repeated the procedure 15 times, each time with a new, dry set of marbles. We found that porosity of the sample is φ = 0.41 ± 0.008 what agrees with literature data [11,40].

Theoretical distribution of the particles forming the granular bed
The YADE code used in further investigations allows introducing data from the above sieve analysis to generate particle clouds with non-uniform diameters. Due to inaccessibility of such data, to use this feature we need to convert the known distribution function to a discrete form with specific number of fractions. Such approach allows generating particle clouds with any distribution. The procedure of generating the theoretical discrete distribution from a Gaussian distribution is as follows. The range of ± 4σ is divided into n f (i.e. number of fractions) equal subranges. Each i-th (where i = 1, 2, . . . , n f ) subrange has an attributed diameter: Then, if N is the total number of spheres in the bed under creation, then the number of spheres n i for the i-th bin is calculated as: Explanation for the above procedure is as follows. The sum of calculated integrals seen in the formula (5) obtained for each bin is equal to one. Multiplication of these values by the number of spheres that are supposed to form the granular bed as a whole, gives us the contribution of each fraction with an attributed diameter.
Note, that the range ± 4σ in case of normal distribution means that 99.994% of the population falls inside that range; in other words, 15,787 must be the minimum total population to expect the first case of value falling outside that range [16].
The result of the procedure described in this section is a cumulative curve of particle distribution given in a discrete form.

Discrete Elements Method (Radius Expansion Method)
To generate a set of virtual beds for which the geometric and hydraulic tortuosity may be next calculated, Discrete Element Method (DEM) [12] implemented in the YADE code [25,41] was used. This method allows modelling the dynamics of physical systems consisting of any number of separate bodies (e.g. particles). During the calculations, two types of equations are solved in a loop. Firstly, the displacement of every body in the system is calculated by the application of Newton's laws, and then, the interaction forces between adjacent bodies are calculated. An important part of the algorithm is tracking contacts between objects. The main equations of the linear and the angular motion may be written as follows [12]: and where A virtual granular bed through DEM simulation may be obtained in three ways: (1) through generation of a particle cloud with proper diameters and size distribution, which is next compressed by moving walls until the target porosity is reached; (2) through generation of a particle cloud with proper size distribution and reduced diameters, that next are (3) through adding particles one by one over the bed, which next fall down by gravity and form the bed. Independent on the approach, all interaction forces between neighbouring particles are calculated during a simulation and particles may move dependent on the local force equilibrium. In our study, the second approach (known in the literature as the Radius Expansion Method, REM) due to the much higher speed of operation was applied (Fig. 2). Some drawback of this approach is that the volume in which the particle cloud grows has to be very precisely defined. In other case, the final porosity or the final average particle diameter would be wrong. This method produces homogeneous beds with slightly varying diameters of particles (due to the differences in the local equilibrium states). We did not take the third method into account due to the conviction that the approach based on particle clouds gives more homogeneous beds. We refer here e.g. to the paper [55] in which this approach was applied in investigations of stratification effects in granular beds.
To calculate the bed volume, so-called characteristic length l c is introduced. This variable is dependent on the number of particles, their size distribution and the whole sample porosity: where and V s is the total volume of the solid part of the porous bed (particles) [m 3 ], V p is the volume of the pore part of the porous bed [m 3 ], V is the total volume of the bed [m 3 ], n j is the number of particles in the i-th fraction [−] and d j the diameter of the i-th particle [m].

Lattice Boltzmann Method
To compute the hydraulic tortuosity with Eq. (3) we need the complete velocity field in pore space [33]. We use Lattice Boltzmann Method (LBM) [7] and the Palabos numerical code [37] to reach this aim. The solver is based on kinetic theory of gases and uses the velocity distribution function f (x, t) in order to describe the state of the system. We solve the transport equation in the form of: where f i (x, t) is the discrete version of the velocity distribution function, c i is one of the lattice vectors, t is simulation time step and Ω i is the collision term. We used the single relaxation time model for collision term (BGK) [47]. The resulting distribution function is time independent. Geometry obtained in a DEM simulation was converted to a grid of nodes needed for a LBM simulation. Some issues related to this step are further discussed.

Path Tracking Method
To calculate the gometric toruosity in virtual beds created by the use of Discrete Element Method, so-called Path Tracking Method (PTM) [43] was applied. The algorithm is suitable for sphere-packed media. In PTM, the key role is played by tetrahedral structures formed by four adjacent particles. Such structures are used to approximate the pore space between spheres. The algorithm finds the shortest path in a given direction (usually chosen as the main direction of the transport in the media) starting from the bottom surface at various locations. The tetrahedrons are constructed one by one and the path length is calculated as a sum of distances between base triangles of all consecutive tetrahedrons. Its ratio to the sample height gives the tortuosity T g [14,43,45]. The algorithm for calculating path length L p may be summarized as follows (Fig. 3): -arbitrarily choose the Initial Starting Point (ISP) at the bottom of the bed; -find three particles (P 1 , P 2 and P 3 ) which are the nearest to the ISP and which form a triangle in the space; -calculate the coordinates of the gravity centre (GC) of the surface (in the triangle plane) through which fluid flows; -move the Initial Starting Point to the Final Starting Point (FSP); in this way the first path section is perpendicular to the bottom surface of the bed; -calculate the normal to the triangle; -estimate coordinates of the so-called Ideal Location (IL), in which the next sphere (P 4 ) surrounding the path should be located; -move the Ideal Location closer to the triangle plane (due to the fact, that particles forming the triangle basis may be separated from one another and in such case the fourth particle is located closer to the triangle plane); -find the particle nearest to the Ideal Location; this is the Real Location (RL) of the 4th particle forming tetrahedron in the space; -remove the lowest sphere from tetrahedron 1-2-3-4 to obtain the base triangle for the next tetrahedron; -calculate the local length of the path inside the current tetrahedron (a local path is the line connecting centroids of base triangles of two neighbouring tetrahedrons); -continue the calculations, until reaching the top surface of the bed.
The algorithm stops when the distance between the gravity centre and the Ideal Location is smaller than the distance between the top surface of the bed and the Z coordinate of the Ideal Location. The last section is perpendicular to the top surface of the bed. The final path length L p is calculated as the sum of all local path lengths inside tetrahedrons and two sections lying at the beginning, as well as at the end of the analyzed pore channel.
The Path Tracking Method was implemented in the PathFinder code [43]. It is a freely available numerical code destined for the analysis of spatial structure of granular beds of a cylindrical or cubical shape, when data about locations and diameters of all particles in the bed are known. Currently two ways of obtaining the input data for the PathFinder code are possible: the indirect (with the use of simulation methods, e.g. DEM) and the direct (with the use of experimental methods, mainly CT/IA). The parameter set calculated by the PathFinder code contains: the bed volume (V ), the volume of the porous part (V p ), the volume of the solid body (V s ), the porosity (φ), the specific surface of the solid body (S 0 ) (with the use of two different definitions), the ratio of actual path length inside channels (L p ), the thickness of the porous medium (L 0 ) and the geometric tortuosity (T g ). Tortuosity may be calculated locally in a chosen ISP or for the whole bed in form of a scalar field [42]. Besides the geometric tortuosity, all other parameters are calculated based on analytical formulas. In the literature, alternative implementations of the PTM may be found [53,54].

Virtual granular beds
All virtual granular beds were created with the use of the Radius Expansion Method described in the Sect. 2.3. Every particle set has a cuboid shape with the size equal to 2l c × l c × l c . The bed size is doubled in the direction in which the geometric and the hydraulic tortuosity is calculated. The procedure of creating a virtual bed was integrated with scripts in Python programing language, needed for running the YADE code (YADE is in fact a set of libraries which has to be called from an other program -here written in Python), and contains following steps: reading the discrete cumulative curve of the particle distribution; calculating the characteristic lenght with the use of Eq. (8); performing DEM simulations; saving the results; converting the data to the other formats (for the Palabos and PathFinder codes).
The generation process was performed for 12 cases, in which the number of fractions n f was being changed. For every fraction, the bed was prepared three times with the use of different values of the seed in the random number generator. In such a way, three independent data sets for every n f were obtained. The number of particles was set to 5000, as this is the number big enough for the reliable comparison of methods of tortuosity determination, allowing performing the analysis in acceptable computation time and not exceeding 15,787, which is the limit for the use of range in the process of generating the theoretical distribution of particles forming the granular bed (as described in the Sect. 2.2).
For the purpose of performing all DEM simulations, the default so-called triaxial test was applied [41]. All settings of the triaxial stress controler were default, besides the final and the maximum multiplier factor of the REM method, which were set to 1.0004 and 1.004, respectively. To obtain a homogeneous bed, the gravity forces were turned off. In order to achieve the final porosity, the friction angle was decreased during a simulation. The starting value of this angle was set to 0.5 [rad]. The time of a typical simulation was about few minutes (Intel Inter Core 3.4 GHz processor, 32 GB RAM and Ubuntu 14.04 OS). An example of a starting and a final particle cloud recorded during a DEM simulation is visible in Fig. 4.
It can be seen ( Table 1) that in every calculation, the obtained bed porosity is slightly different, even for the cases in which the number of particle fractions was constant. This is a characteristic feature of the used method. Note that the characteristic length is a little different in every case.
The porosity in every case is a little less than the target (experimental) porosity due to the condition of the calculations termination-the loop is finished when the current porosity is less than the target porosity. The relative errors (δφ) between the obtained porosities and the experimental value are in all cases less than 1%. For this reason we can assume that the generated virtual beds represent the real bed sample very well.
Other important issue is to check the conformity of diameters in the real and the virtual bed. In Fig. 5 three examples of the diameter distribution (for n f = 3, 13, 25 and data set 1) of particles included in the virtual bed, as well as the normal distribution, are shown. It can be seen that diameters in the virtual bed not coincide exactly with the theoretical distribution prescribed based on measurements of the real particles. It is the result of the way of generation of the virtual bed. The mean diameters of the spheres in generated beds are as follows: 6.0717 (for n f = 3), 6.0774 (for n f =13) and 6.0692 (for n f = 25). As the average diameter of theoretical distribution is equal to 6.0720, it can be stated that the differences are negligible (0.0043%, 0.0884% and 0.0463%, respectively). The Pearson's linear correlation coefficients [39] were also calculated between histograms for theoretical distribution and histograms for distributions of obtained virtual beds. The results (99.99% for n f = 3, 97.81% for n f =13 and 98.08% for n f = 25) are satisfactory here as well. The same results were obtained for others data sets.

Hydraulic tortuosity (LBM simulations)
To perform calculations in the Palabos code (see Sect. 2.4), the bed geometry obtained in a DEM simulation was transferred to the form required by the LBM. A regular, cubic lattice of nodes was used to approximate geometry and porosity of the discretized sample. The size of each grid was equal to 2L × L × L. Due to the high demand for computing power only three cases of n f were analysed.
To find the minimal system size we discretize geometry of each sample and set each node of the lattice as solid if it lays entirely or at least by 50% inside any of the beads. Remaining nodes are left as fluid. The ratio between a number of fluid nodes to all lattice nodes gives porosity. We test various values of L ∈ 75, 100, . . . , 300 [lu] (lattice units) and find the minimal size L = 200 [lu] at which porosity does not change with an increasing grid size (Fig. 6).
In the next step, we measured hydraulic tortuosity using the same samples. First, we simulated the fluid flow at the pore-level. We set periodic boundaries in the flow direction. The no-slip conditions were used at horizontal walls. Also, the no-slip conditions were used for solid nodes. The relaxation time was τ = 1 resulting in the kinematic viscosity ν = 1/6 in lattice units [38]. The fluid was driven by the external gravity and the creeping flow was kept (the Reynolds number Re= dq/ν < 1, where q is a Darcy flux, d is the mean grain diameter). An exemplary result of the simulation is given in Fig. 7. The data obtained for L = 200 are collected in Table 2. Hydraulic tortuosity was calulated using Eq. 3. To check the finite-size effects we gradually increased the lattice size from L = 75 to L = 300. We found (Fig. 8) that tortuosity grows monotonically following an exponential function. Thus, in order to estimate the value of T h for infinite lattice, we fit an exponential function of the form where T (∞) = 1.24 and fit parameters a and b are equal to 0.13 and − 0.014, respectively. Details of the procedure are similar to those used i.e. in [32]. Please note that T h ≡ T (∞) from the above mentioned fitting procedure.

Geometric tortuosity (PTM calculations)
The data obtained in DEM simulations was converted to format readable by PathFinder code (see Sect. 2.5), which was next used to calculate the porosity and the geometric tortuosity. In all calculations default settings were used [43].
The results of calculations are collected in Table 3. It can be seen (Fig. 9) that the porosity in every case is a little higher than the experimental porosity. The relative errors are usually slightly higher than these obtained in the direct calculations in the YADE code. However, we assume that these values are still very well representative. This is all the more justified that the differences in the porosity between particular ways of its calculation are less than the measurement error (φ ex p ± Δφ ex p ). The geometric tortuosity was calculated 25 times for every virtual bed in the direction of the longest side of the calculation domain (in Table 2 only the average values are shown). Such a large number of values results from the fact that the path lengths inside channels were calculated in different zones of the virtual bed (Fig. 10) [43]. In this way, an average value was obtained, which better represents the tortuosity in this bed. In Fig. 11 the average values of geometric tortuosity for all three data sets and different particle fractions are shown. The average value of the tortuosity obtained with the use of PTM for all data sets and all values of n f is equal to 1.2147. In this figure, the results reported by Wang (see Sect. 1) are also shown. The relative error between the above mentioned average value and the results shown by Wang is equal to 1.43 %. In such a way, the usefullnes of the PTM is validated not only by LBM simulations, but by an independent results as well.
No systematic changes with n f were found, but, for n f > 10 values of T g are found to be more scattered, higher and encumbered with larger statistical errors of the mean.

Direct data comparison
In our simulations we show that the porosity in all used numerical models is very close to the experimental value. All the realative errors were close to 1% or (usually) even less than 1%. The particle distribution obtained by the use Table 2 Results of calculations in the Palabos code for L = 200 (data averaged over three independent samples for chosen n f ) Influence of the numerical grid resolution on the hydraulic tortuosity of REM follows the theoretical distributions at a satisfactory level. We noted that the higher values of the n f give more representative virtual beds due to the better reconstruction of the Gaussian distribution. Finally we can assume that if the porosity value is correct, than the other, more complicated and porosity dependent, physical quantities may be also determined from the same virtual model (e.g. tortuosity). We observed that average values of hydraulic and geometric tortuosity are very similar and slightly higher than the results shown by Wang. Probably the reason is that Wang used particles with the same diameter. It suggests that tortuosity may increase with the increase of the standard deviation of the particle distribution. However, to confirm this, more investigations are needed.
The main difference between the both used methodologies is the scattering of the obtained tortuosity values (Fig. 12). In the hydraulic approach, the tortuosity value is averaged over all nodes (available for the fluid) of the lattice grid. The   Fig. 9 Comparison of the porosity value obtained in the all used methods local values of the tortuosity are not known. For this reason the average value is not very sensitive to the arrangement of the particles in the bed (i.e. the used data set). We suppose that the difference between average tortuosities obtained for different data sets would decrease with increase in the number of particles in the bed. This feature may be treated as an advantage (we obtain a very similar value every time) or as an drawback (we have no information about the possible local fluctuations of the tortuosity). In this context, the geometrical approach gives more information (e.g. in paper [42] twodimensional tortuosity fields in the direction perpendicular to the main flow direction were calculated). We obtain not only one average value, but also the data on the range of the geometric tortuosity.
It should be underlined that the range, as well as the average value of the geometric tortuosity, is in some way dependent on the number of the Initial Starting Points used in calculations (see Sect. 2.5). In the presented study, always 25 points were used. But, if we change the number of the   Fig. 11 Relationships between the number of fractions and the tortuosity

Indirect data comparison
The observed relative differences between geometric and hydraulic tortuosity in granular beds may be expressed in terms of the anisotropy of an inclined capillary model of the flow [6]. In that model particles paths follow straight paths inclined by some angle α (varying from path to path) with respect to the main axis of the flow. The tortuosity of each channel is equal to where α is the inclination angle. To describe the whole medium, L e f f may be taken as the length of an average flow axis. Using the above relation and known values of T g and T h the inclination angles for geometric and hydraulic tortuosity are evaluated α g = 34 • and α h = 36 • , respectively. From this, the absolute difference 2 • and the relative difference Δα = 3 • is calculated. Although this analysis is not perfect and suitable for simplified model of inclined channel flow, it gives some rough idea how big error would be introduced by taking T g instead of T h for the anisotropy angle estimations. In the real porous medium, due to a non-straight, wiggly paths of particles, one may expect that this difference will be decreased. Similarly, if we consider the empirical law for permeability (k) as a function of geometric parameters and tortuosity in which k ∝ 1/T 2 [23,51], the maximum deviation Here k g is based on geometric and k h is based on the hydraulic estimation of tortuosity. The same 6% deviation appears between effective diffusion coefficients in porous media related to diffusion in fluids with D e f f = DT −2 [8].
Our results may be also used to analyze sphere-packed porous media in terms of the Archie law. This is an empirical relation for tortuosity as function of porosity: where φ is porosity and b is the free parameter related to a cementation factor in porous media [5,10]. Using Eq. (14), for porosity φ = 0.41 (Sect. 2.1) and tortuosity 1.24 (Sect. 3.2) and 1.21 (Sect. 3.3) we get b ≈ 0.24 and b ≈ 0.21. This is close to b = 0.25 calulated previously by numerical fitting of the Eq. 12 to tortuosity of the random, overlapping spheres three-dimensional porous medium [33]. This is an unexpected result which suggests that the cementation factor in the Archie law is independent of the way spheres are packed and arranged. We suggest that this is related to the shape of particles packed in the bed. This finding deserves further study and analysis but it is out of the scope of this work.

Possibility of investigations of the bed anisotropy
In the paper, we discuss results of tortuosity calculations in only one space direction. However, both used methods may be used to determine the tortuosity in others directions. Having such data, any three-dimmensional model destined to analysis porous media may be used. The so-called Porous Media Model available in the popular ANSYS Fluent code may serve here as an example. The linear (Darcy) term in this model has a form [1]: . Now, if we assume that the term D is described e.g. by the Kozeny-Carman law (1), we can predict pressure distribution in a three-dimensional computational space. Note that this law does not define, which kind of the tortuosity should be used. We are aware that the anisotropy effects in granular beds consisting of spherical or quasi spherical particles are usually not very significant. However, in beds with high standard deviation, the particle diameters may be very different and the particles may tend to aggregation or form layers due to sizes. The way of generating such porous beds (real or virtual) may also affect in the level of the anisotropy. Except that, in a general case, the particle distribution may not follow Gaussian distributions and may by a mixture of many fractions with different distributions. The PTM may be applied for porous beds of each above kind.

Time and memory consumption
We observed that computation time used to calculate the geometrical tortuosity was signicantly less than the time needed for LBM solver to simulate fluid flow required for T h computation. For example, to determine single path in geometric methods using the earlier mentioned computer less than 6 [s] was necessary in average. In LBM, computation of fluid flow at reasonable accuracy required for T h measurement varied with mesh size L and it was less than 20 [min] for L = 75, more than 1.5 [h] for L = 200 and more than 5 [h] for L = 300. Simulations in Palabos were performed in parallel mode with the use of 4 processors. Similarly, the memory consumption of LBM was large. It was 3GB for the D3Q27 model at L = 200 and 11GB for L = 300. For the same sample, PathFinder stores only grain positions and diameters, which required less than 0.2 MB for 5000 beds. Even using the sparse data structures in Palabos, that could drop the memory down to 45% [29], the memory required by LBM is orders of magnitude larger than required by geometrical algorithms.
By taking into consideration lower time and memory resources required by geometric methods, we suggest that geometric methods will allow to estimate tortuosity in much larger, granular systems. This may also simplify measurement of other important parameters (i.e. permeability) in large scale samples.

Conclusions
We found that geometric tortuosity in granular beds with slightly varying diameter is T g ≈ 1.21, which is close to the hydraulic counterpart T h ≈ 1.24. This value is also close to the hydraulic tortuosity given by Wang, who used a similar methodology of obtaining this quantity. We suggest that the applied geometrical method may serve as an attractive alternative to hydraulic tortuosity in large granular systems. The advantages of the Path Tracking Method are as follows: a very short time of calculations; the possibility of use in granular beds consisting of a large number of particles; the possibility of calculating tortuosities of individual geometric path lines in one virtual bed; no need of geometry conversion; no need of model calibration. In turn, in the LBM approach the demand for computing power is large, what limits its use.
The main drawback of the PTM is that it is not related directly to the transport equations and that this method is appropriate only for beds consisting of spherical or quasi-spherical particles. However, the area in which such kinds of particles occurs is still relatively large.