A mathematical model of tortuosity in soil considering particle arrangement

Tortuosity is an important parameter for studying the permeability of soil. Existing studies of soil tortuosity are usually of empirical nature and attempt to relate tortuosity to soil porosity alone. By assuming a laminar flow through the pores of two‐dimensional square solid particles, we present a mathematical model for calculating soil tortuosity under different particle arrangements. The effect of the randomness of the particle arrangement on the tortuosity is evaluated, which generates the variation range of the tortuosity. The proposed model provides the upper and lower bounds of the tortuosity, while existing empirical models tend to fall within these bounds. The consistency between the proposed model and the numerical calculation provides a validity for the proposed model.


INTRODUCTION
The soil hydraulic conductivity is of great importance in the study of the percolation process in porous media (Bear, 1988;He, Teng, Sheng, & Sheng, 2017;Sheng, Zhang, Niu, & Cheng, 2014;Teng, Kou, Zhang, & Sheng, 2019a), consolidation and settlement of soils and foundations (Ren et al., 2016;Rodriguez, Giacomelli, & Vazquez, 2004;Zhang, Sheng, Zhao, Niu, & He, 2015), migration of pollutants from waste disposal facilities (Malusis, Shackelford, & Olsen, 2003;Yanful & Choo, 1997), and other geotechnical engineering problems. Theoretical and empirical equations relate hydraulic conductivity to soil properties like soil texture and structure, pore geometry, and chemical properties of the pore fluid (Costa, 2006;Ren et al., 2016;Teng, Shan, He, Zhang, & Sheng, 2019b;Teng, Zhang, Zhang, Zhao, & Sheng, 2019c). Brooks and Corey (1964, pp. 352-366, Abbreviations: LA, lower limit arrangement; UA, upper limit arrangement. This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. (1) where v is the averaged flow velocity (m s −1 ), n is the porosity (unitless), R is the averaged hydraulic radius of the pores (m), c is a shape factor of the pores (unitless), μ is the dynamic viscosity of the fluid (Pa s), ∆P is the pressure difference between the two sections (Pa), L is the height or length of the sample (m), and Γ is the tortuosity of the soil (unitless) which is the ratio of the actual flow path of the fluid L t (m) to parameter L (Yuan, 2008). By combining Equation 1 with Darcy's law, we obtain the expression for hydraulic conductivity, k: Vadose Zone J. 2020;19:e20004. wileyonlinelibrary.com/journal/vzj2 1 of 12 https://doi.org/10. 1002/vzj2.20004 where k is hydraulic conductivity (m s −1 ), ρ is the density of the fluid (kg m −3 ), and g is the acceleration due to gravity (m s −2 ).
As described in Equation 2, the tortuosity plays a crucial role in describing k. The concept of tortuosity was introduced by Carman (1937Carman ( , 1939Carman ( , 1956 to match the hydraulic conductivity computed based on a bundle of capillary tubes to experimental data (Dullien, 1979;Zhang et al., 2016). The hypothesis derived from Carman (1956) is that the pores in the soil are regarded as curved "capillary tubes". Fluid can flow in these tubes. These "tubes" in the soil have three key factors: length of tubes, area of tubes, and number of tubes. Just as the shape factor c describes the shape of the cross-section of the "tubes," the effect of the tortuosity is to describe the actual length of the "tubes." The concept of tortuosity seems straightforward and is meant to approximate the actual length of the paths that fluid molecules travel through the pore pace (Ghanbarian, Hunt, & Ewing, 2013a). In reality, however, tortuosity is often used as an adjustable parameter, i.e., a "fudge factor," which implies the lack of a clear understanding of tortuosity.
It has been recognized that the tortuosity is positively related to the porosity of the soil. Various empirical expressions have been presented to describe this relation in the literature, as shown in Table 1. Comiti and Renaud (1989) presented an empirical logarithmic function to relate tortuos-T A B L E 1 Empirical expressions of tortuosity (Γ)

No. Equation Source
1 Γ = 1 − ln( ) Comiti and Renaud (1989) 2 Γ = 1 − 0.8(1 − ) Koponen et al. (1996) 3 Koponen et al.  Yu and Li (2004) Note: P, a, and m are fitting coefficients, with a usually set to 0.65 and m to 0.19, n is the porosity, n c is the threshold porosity and equals 0.33, D p is the solute diffusion coefficient in soil (cm 2 ) as defined by Fick's law of diffusion, D 0 is the solute diffusion coefficient in water (cm 2 ), f is the volumetric fluid-phase content, ξ is the sphericity (or roundness) factor, θ is water content, θ t is threshold water content, υ is a scaling exponent, D x is the fractal dimensionality, L s is the straightline length across the medium, and C is a numerical factor.

Core Ideas
• The effect of particle arrangement on the tortuosity is evaluated. • A new mathematical model for computing tortuosity accounts for different particle arrangements. • The model provides upper and lower bounds of tortuosity; existing models fall within these bounds.
ity Γ to porosity n. Comparing their results using this function with the measured results of particle bed fluid experiments, they indicated that a fitting coefficient P was equal to 0.41 for spherical particles and 0.63 for cubic particles. Mauret and Renaud (1997) noted that P should be 0.49 in a capillary model of a high-porosity bed composed of spheres and fibers. Koponen, Kataja, and Timonen (1996) simulated the flow of incompressible Newtonian fluid through freely arranged square particles by means of the automatic lattice gas method and derived a linear relationship between Γ and n. Koponen, Katya, and Timonen (1997) further considered the presence of the permeability threshold and proposed a new power function. Moldrup, Olesen, Komatsu, Schjonning, and Rolston (2001) developed an empirical formula for tortuosity based on the solute diffusion process. Some more theoretical models of tortuosity have been developed based on geometric properties (e.g., particle size, shape, and arrangement) and topological properties (e.g., the dimensionality and connectivity of the network) of porous media (Ghanbarian et al., 2013a). For example, Yu and Li (2004) proposed a tortuosity model by assuming two-dimensional square particles in an equilateral-triangle arrangement. The formula has no empirical parameters and can explain the relationship between tortuosity and porosity, as shown in Table 1. Lanfrey, Kuzeljevic, and Dudukovic (2010) assumed that each flow path was represented by a sinuous tube with constant cross-sectional area and perimeter and developed a theoretical model for the tortuosity of a fixed bed. The particles were positioned in the form of an equilateral triangle and square unit cell arrangements. Li and Yu (2011) obtained a tortuosity model based on the hierarchical structure of a deterministic Sierpinski carpet, which belongs to a pore fractal model that includes particles of different sizes and pores of the same size (Rieu & Perrier, 1998). Based on percolation theory and the finite-size scaling approach, Ghanbarian, Hunt, Sahimi, Ewing, and Skinner (2013) proposed an expression of tortuosity that is applicable to saturated and unsaturated porous media. In fact, the actual arrangement of soil particles is very complicated. Even if the macroscopic porosity of the soil is the same, the arrangement of particles can differ considerably, which causes variation in the flow path and tortuosity of the soil (Ghanbarian et al., 2013a). Some studies in the literature found that tortuosity is really related to particle arrangement (Ho & Strieder, 1981;Tsai & Strieder, 1986;Comiti & Renaud, 1989). However, the relationship between tortuosity and particle arrangement remains unknown. The simple assumption about particle arrangement should be overcome to allow a new understanding of tortuosity in soil.
The objective of this study was to identify the effect of particle arrangement on tortuosity and to establish new geometric tortuosity models considering different particle arrangements. Furthermore, the macroscopic and microscopic characteristics of tortuosity are explained by comparing with existing models. The COMSOL Multiphysics package was used to simulate the soil and its pore fluid, and the results of the numerical simulation were compared with the present model for verification.

MATHEMATICAL MODEL OF TORTUOSIT Y
To derive a mathematical model of tortuosity, we made the following assumptions: (a) soil particles are represented by two-dimensional squares of identical size, which is consistent with the assumption adopted in the literature, for example by Yu and Li (2004); and (b) the pore fluid in the soil is a Newtonian fluid, and its flow in the soil is assumed to be laminar. The soil particles are evenly distributed as shown in Figure 1a, where A is the length of the soil particles and B and C are the distances between two soil particles in the vertical and horizontal flow directions, respectively. An anisotropic parameter m is defined here to describe the ratio between B and C, i.e., m = B/C. The retarding parameter θ is the horizontal angle between two particles, which ranges from 0 to arctan(B/2C). The special case when θ equals 0 is defined as the lower limit arrangement (LA), as shown in Figure 1b. The normal case when θ equals arctan(B/2C) is defined as the upper limit arrangement (UA), as shown in Figure 1c. The parameter θ indicates the blocking effect of soil particles to the liquid water flow. Therefore, the anisotropic parameter m and retarding parameter θ determine the relative position of the particles, i.e., the arrangement of the soil particles.
The soil porosity n can be determined by Equation 3a, which can be further transformed into Equation 3b: where V t is the total area of the selected region, which is the parallelogram area in Figure 1a, and V s is the area of the particles in the selected region. Since the anisotropic parameter m is the ratio of B to C, we can obtain the relationship between A and C: In this model, it is assumed that A must be smaller than C and B, which means that The assumed laminar flow in this study can be described as the flow path characterized by many parallel lines, as shown in Figure 2a. To quickly compute the flow path, three regions can be defined according to shape, as shown in Figure 2b. The flow in Regions I and III represent the liquid flows around particles, and the flow in Region II represents the liquid flows in the pore. The flow path in Regions I and III is A/2. In Region II, some topical flow path is chosen to compute. The longest flow path is the flow path between a 1 and d 1 , and the shortest path is the flow path between a 2 and d 2 . The length of the flow path between a 1 and d 1 ranges from 1 1 to 1 1 1 1 , and their average is defined as L max : The shortest path, L min , can be expressed by Considering that the flow path becomes smoother when the porosity increases, Equations 5a and 5b can be transformed into According to Equations 6a and 6b, the averaged flow path in Region II is Vadose Zone Journal F I G U R E 1 (a) Schematic diagram of the particle arrangement, (b) the lower limit arrangement (LA), and (c) the upper limit arrangement (UA). Here, A is the length of the soil particle and B and C are the distance between the two soil particles in the vertical and horizontal flow directions, respectively. The parameter θ is the horizontal angle between two particles F I G U R E 2 (a) The assumed laminar flow paths, where A is the length of the soil particles, B and C are the distances between two soil particles in the vertical and horizontal flow directions, respectively, the anisotropic parameter m is the ratio between B and C, and the retarding parameter θ is the horizontal angle between two particles; and (b) the division into regions to compute the flow path. The flow in the Regions I and III represent the liquid flows around particles, and the flow in Region II represents the liquid flow in the pores The particles will overlap in the actual situation. If the particles overlap, the flow paths in Regions I and III can be ignored when we compute the total flow path. Only the flow paths in Region II are considered. In this case, the tortuosity of the soil is Γ 1 , which can be expressed as the ratio of the flow path in Region II to length C − A: If we neglect the particle overlap, the soil tortuosity Γ 2 can be computed as the sum of the flow paths in Regions I, II, and III and length C: Substituting Equation 4 into Equations 8a and 8b leads to the final expression of Γ: Taking an average of Γ 2 and Γ 1 , we obtain the expression for tortuosity Γ as Equation 10 considers the effect of both soil porosity and particle arrangement on the tortuosity. The soil particle arrangement is controlled by the anisotropic parameter m and angle θ. The particle arrangement is considered less to compute tortuosity Γ in the literature. Through analysis, the parti-cle arrangement affects the flow path of the fluid. Equation 10 shows that the tortuosity has an upper limit and a lower limit due to the effect of the particle arrangement, which correspond to UA and LA, respectively. Equation 10 provides a new method to compute the tortuosity based on the porosity and particle arrangement, which is more universal and applicable.

VERIFYING THE PROPOSED TORTUOSIT Y MODEL
To validate the proposed model, a numerical approach was adopted to simulate the hydraulic conductivity; then, the tortuosity was calculated according to the relationship between hydraulic conductivity and tortuosity. In this study, the COM-SOL Multiphysics program package was used to simulate water flow through the square particles of soil as the laminar flow of a Newtonian liquid. As shown in Figure 3, square particles represent the soil particles, the blue part is the fluid Inlet/outlet length, mm 8 The length of B, mm 2

Reynolds number <5
Note: B is the distance between two soil particles in the vertical direction; a Reynolds number of less than 5 is to guarantee the applicability of Darcy's law.
portion, and the white arrows indicate the direction of flow. The boundary of the wall is a symmetrical boundary, and the system is designed to be periodic, which means the system is homogenized. The fluid, with a constant initial velocity, continuously flows in from the left inlet, flows through the particles, and finally flows out through the right outlet. Some probes are placed in the left and right end to test the flow rate and water head difference. The model considers two types of limit arrangements: LA and UA. The values of A, C, and θ can be changed to simulate different particle arrangements. The settings and the equations used in the numerical simulation are shown in Tables 2 and 3, respectively. The COMSOL Multiphysics program package solves the equation by the finite element method. The solver used in this study was the Multifrontal Massively Parallel Sparse Direct Solver (MUMPS), which is a direct solver based on LU decomposition. The COMSOL Multiphysics program package uses physics-controlled meshing sequences and free triangular mesh to generate a numerical grid. The mesh quality and mesh resolution are shown in Figure 4.
Relative tolerance was adopted as the stop criterion to solve the time-dependent questions. If the relative residual in the solution is lower than the relative tolerance, the solution will continue, and the result will converge.
It is noted that the numerical simulation model is scale dependent, and the models should satisfy the Reynolds similarity principle. The Reynolds number is defined as where d is the characteristic linear dimension (m). In this model, the parameter v is the inlet flow velocity, and d is equal to the side length of particle A. The same Reynolds number can ensure the similarity among different cases. Therefore, it is important to select an appropriate Reynolds number for the computation. Reynolds numbers were set as 2, 0.2, and 0.02. The calculation parameters are shown in Table 4. The computed results are shown in Figure 5. It can be seen that the same hydraulic conductivity was generated corresponding to the same porosity when te Reynolds number was 0.2 or 0.02. When the Reynolds number was 2, the hydraulic conductivity corresponding to a high porosity was smaller than the others. Therefore, the Reynolds number was assigned to be 0.2 in the following numerical simulation. According to Darcy's law, the hydraulic conductivity of the soil can be expressed as = ρ s Δ Δℎ (12) where Q is the flow rate (m 3 s −1 ), ρg is the gravity of the fluid (N m −3 ), A s is the soil area (m 2 ), ∆l is the soil height (m), and ∆h is the head difference (N m −2 ). The parameters A s and ∆l are shown in Figure 3. By substituting the parameters, such as the flow rate, water head difference, and so on, into Equation 12, the hydraulic conductivity for different particle arrangements can be obtained.
Considering that the numerical simulation model is a planar, two-dimensional, parallel plate model of laminar flow, Equation 2 can be used to express the relationship between the hydraulic conductivity and the tortuosity: where n' is the ratio of the cross-sectional pore area to the cross-sectional area (unitless), and b is the interparticle void size (m), as shown in Figure 3. Parameters n', b, and Γ are affected by the porosity and particle arrangement. Carman

Condition Equation Value
Governing equation Note: ρ is the density of the fluid, u is the velocity matrix of flow, t is time, p is pressure, F is the mass force, μ is the dynamic viscosity, I is a unit tensor, n is a unit vector, V 0 is the initial velocity, Re is the Reynolds number, d is the characteristic linear dimension, and p 0 is the pressure at the outlet. Note: A is the length of soil particles, Re is the Reynolds number, V 0 is the initial velocity, μ is the dynamic viscosity, ρ is the fluid density, d is the characteristic linear dimension, m is the anisotropic parameter, and θ is the retarding parameter. (1937,1939,1956) and Duda, Koza, and Matyka (2011) found that the shape factor is independent of the porosity. Therefore, the same particle arrangement is assumed to have the same shape factor. In this study, we regarded the shape factor c as the value when the porosity equals 0.5.

Comparison between the proposed model and the numerical simulation
In numerical simulation, the change in porosity n is controlled by changing the particle length A and the change in m is controlled by changing the value of C. The inputs for the models are presented in Tables 5 and 6.

F I G U R E 5 Comparison of the hydraulic conductivity for different
Reynolds numbers (Re) Because in the model there is no particle overlap, the theoretical result of the tortuosity is determined by substituting the corresponding parameters of the numerical simulation into Equation 9b. The simulated solution and theoretical solution were compared to verify the rationality of the model, as shown in Figure 6. The tortuosity obtained by the theoretical solution decreases with increasing porosity, and the numerical Note: A is the length of soil particles, Re is the Reynolds number, V 0 is the initial velocity, μ is the dynamic viscosity, ρ is the fluid density, d is the characteristic linear dimension, m is the anisotropic parameter, and c is a pore shape factor.

7.95
Note: A is the length of soil particles, Re is the Reynolds number, V 0 is the initial velocity, μ is the dynamic viscosity, ρ is the fluid density, d is the characteristic linear dimension, θ is the retarding parameter, and c is a pore shape factor.
simulation results have the same trend. The results indicate that the particle arrangement affects the tortuosity, and the effect decreases with increasing porosity. The results of the theoretical calculation and numerical simulation are highly consistent with each other when the porosity ranges from 0.3 to 0.9, which proves that the proposed model has a relatively high accuracy. It was also found that there is a gap between the theoretical result and numerical simulation when the porosity is too high or too low. This is largely due to the difference between the actual streamline and that assumed in the mathematical model. The gap was mentioned by Matyka, Khalili, and Koza (2008) as well. They concluded that the error can be caused by a flaw in the numerical simulation method, for example, the difficulty of reaching the stationary solution and the existence of discontinuities in stream line. In this study, the method of statistical streamline was not used to calculate tortuosity, which may avoid part of the errors. Another reason for this error is the difference between the actual streamline and that assumed in the mathematical model. Because tortuosity is the ratio of the length of the actual flow path to the sample length, the minor difference in numerical and theoretical results does not change with sample size.

Comparative analysis with the methods in the literature
The predicted model in this study was compared with methods in the literature, including those of Comiti and Renaud (1989), Koponen et al. (1996Koponen et al. ( , 1997, Yu and Li (2004), and Ghanbarian et al. (2013). The computed results of these models are shown in Figure 7. It can be observed that most of the results of the tortuosity model in the literature fall within the region of the proposed model. The proposed model provides the upper and lower limits for the tortuosity value, which offers an insight iinto understanding tortuosity. It indicates that the particle arrangement has a remarkable effect on the tortuosity. The proposed model can provide a variation range for the tortuosity instead of a certain value, which is considered more reasonable. The results in Figure 7 also indicate that the particle arrangement in the UA always results in a larger tortuosity than the LA. In the case of the UA, the flow path will be greatly increased; thus, the fluid is flows with difficulty. When the porosity is greater than 0.5, the tortuosities of the two arrangements converge and finally reach a value of 1.0 when the porosity is close to 1.0. The result also indicates that when the porosity is large, the particle obstruction to the fluid will be weakened, and the fluid can more smoothly pass through the particles. The particle or pore size in the proposed model refers to the relative size of the particles or pores and not the actual size. When the porosity is high, or the particles are smaller than the pores, the particle size hardly affects the tortuosity. In this case, most of the flow paths can be approximated as straight lines, so the tortuosity is nearly 1. When the porosity is low, the particles are large relative to the pores. The particle size greatly contributes to the tortuosity, and the degree of tortuosity sharply increases. We note that the result of Comiti and Renaud (1989) lies below the LA when the porosity becomes lower. This may be caused by the assumption in their study of a rectangular particle, while a square particle was used in this study.

PARAMETRIC ANALYSIS
There are three parameters in the proposed model of the tortuosity: porosity n, anisotropic parameter m, and retarding parameter θ. The void ratio is mainly determined by the size of the particles and the distance between the particles, which are macroscopic quantities. The anisotropic parameter m and retarding parameter θ mainly determine the particle arrangement and can significantly affect the tortuosity, which will be discussed in detail.
The anisotropic parameter m reveals that the tortuosity is associated with the seepage direction. In the above assumption, the fluid flows into the particle from side B. However, the tortuosity of the fluid flowing into the particle from side B is different from that from side C. In Equation 10, changing the flow direction of the fluid is equivalent to changing the size of the anisotropic parameter m. When the particles are arranged at the lower limit (θ = 0), the particles are isotropically arranged, i.e., m = 1. When the particles are arranged in anisotropy, it is assumed that the fluid flows into the particles from side B, and m = 0.8 (i.e., B = 0.8C). If the fluid flows into the particles from side C, m = 1.25 (C = 0.8B).
Substituting the above values of m into Equation 10, we have the calculation result in Figure 8a. When porosity n is small, a greater m value always leads to a greater tortuosity. With increasing n, the tortuosity tends to 1, and the effect of m on the tortuosity becomes unapparent. When n = 0.25, the tortuosity can be 1.5 for m = 0.8 or 4.5 for m = 1.25. This value reflects that the flow from one direction in the seepage is difficult, but it is smooth from the other direction. The results indicate that the change in seepage direction can lead to different tortuosities even for the same soil, which results in different hydraulic conductivities.
The retarding parameter θ also controls the particle arrangement. If parameter m is equal to 1, parameter θ is assigned to be 5, 15, and 25 • , and the relation between tortuosity and porosity is shown in Figure 8b. A greater θ always leads to a greater tortuosity because θ indicates the blocking effect of soil particles to the liquid water flow; a greater θ results in a longer flow length.
It is noted that the goal of this study was to analyze and describe the relation between particle arrangements. The proposed model can take the particle contact, porosity changes, existence of anisotropy, and heterogeneity of the system into account, which can be regarded as a simplified model for soil particles.
Particle contact: the proposed model takes the influence of particle overlap into account.
Porosity variation: porosity is one of the input parameters in this model. Anisotropy: the anisotropic parameter m is used in the proposed model, which can explain the anisotropy of the soil.
Heterogeneity: heterogeneity refers to the uneven arrangement of particles at the same porosity. The limit arrangement can be determined as well as the value range of tortuosity at the same porosity.
However, some limitations still exist in the proposed model, for example: 1. The proposed model is a two-dimensional model, and the particles are square shaped and well arranged, while soil particles in actuality are three dimensional, with irregular shape and random arrangement. The assumption should be satisfied that A is smaller than both C and B. When A is larger than B, it means that all the channels in the direction of fluid inflow are blocked. When A is larger than C, the hypothesis of the flow path differs greatly from the actual situation. Based on this assumption, the threshold porosity can be determined as n threshold = max{1 − m, 1 − 1/m} in this model. These cases should be studied further.

CONCLUSION
Tortuosity is important in understanding the permeability process of soil. However, there has been no generalized formula for calculating the tortuosity due to various complex factors. In this study, a new mathematical model of tortuosity is established by considering the particle arrangement. The main findings are: 1. The proposed model provides the upper and lower bounds of the tortuosity, which correspond to two types of particle arrangement: upper and lower limit arrangements. By comparing the results of the model with the results of analytical or empirical methods in the literature, we show that the proposed model can provide the variation range of the tortuosity, and the results of most models fall in this range. Moreover, the proposed model is consistent with the numerical calculation, which proves the validity of the proposed model. 2. The parameter analysis shows that the particle arrangement significantly affects the tortuosity when the soil is relatively tight (low porosity), but this effect becomes less pronounced when the soil is loose (high porosity).
3. The proposed model is two dimensional, and it should be extended to three dimensions in future study.

NOTATION
average flow path in Region II, m A length of soil particles (See Figure 3), m A s soil area (see Figure 3), m 2 B distance between two soil particles in the vertical flow direction (see Figure 3), m b interparticle void size (see Figure 3), m C distance between two soil particles in the horizontal flow direction (see Figure 3), m c shape factor of pores, unitless d characteristic linear dimension, m g acceleration of gravity, m s −2 k hydraulic conductivity, m s −1 L height or length of sample, m L max longest flow path in Region II, m L min shortest flow path in Region II, m m anisotropic parameter, defined as the ratio between B and C, i.e., m = B/C n porosity, unitless n ′ ratio of the cross-sectional pore area to the crosssectional area, unitless Q flow rate, m 3 s −1 R averaged hydraulic radius of pores, m Re Reynolds number, unitless ∆h head difference, N m −2 ∆l soil height (see Figure 3), m ∆P pressure difference between two sections, Pa v averaged flow velocity, m s −1 Γ tortuosity of the soil, unitless Γ 1 tortuosity considering particle overlap, unitless Γ 2 tortuosity neglecting particle overlap, unitless θ retarding parameter, defined as the horizontal angle between two particles (see Figure 3), unitless μ dynamic viscosity of fluid, Pa s ρ density of fluid, kg m −3