Coarse-grained geological porous media structure modeling using heuristic algorithm and evaluation of porosity, hydraulic conductivity, and pressure drop with experimental results

Knowledge of porous media structure is an essential part of the hydrodynamic investigation of fluid flow in porous media. To study soil behavior (as a granular porous media) and water and contaminant movement in the vadose zone, appropriate estimation of soil water retention curve (SWRC) and soil hydraulic conductivity curve (SHCC) has a pivotal role and is one of the most challenging topics for researchers and engineers in soil and water science. The SWCR can be approximated using an accurate particle size distribution (PSD) function. In this study by applying the random close packing (RCP) method as an encouraging method for predicting and studying particle configuration, an optimal particle size distribution is developed for coarse-grained soils (0.025 mm < PSD < 3.35 mm). The mentioned RCP is generated using a heuristic algorithm with merging applicable equations of soil science. For porous media modeling, MATLAB software is used and the predicted results by the optimal model for the parameters of porosity, pressure drop, and saturated hydraulic conductivity are compared with laboratory measurements. Experimental design is conducted by MINITAB and predicted coarse-grained soil structure by the model is compared with four sifted soils. The results of the sensitivity analysis showed that the porosity obtained from the model is strongly sensitive to the resolution factor and should be chosen with a sufficiently large amount (higher than 250). Results showed good consistency (up to 95%) between predicted porosity and only a 10% difference in pressure drop and permeability with observed measurements.


List of symbols d
Diameter (

Introduction
Soil water retention curve (SWRC) is a powerful and instructive tool to prognosticate soil behavior and water and contaminants movement in the vadose zone (Tong et al. 2012). SWRC describes the soil matric suction (h) as a function of volumetric soil moisture (θ) or vice versa and along with soil hydraulic conductivity curve (SHCC) is utilized to estimate soil hydraulic properties (Braddock et al. 2001). There are several empirical and theoretical methods for predicting SWRC in literature. The ASTM D6836-16 has been developed to laboratory measure SWRC consisting of five standard methods (A-E) at various suction ranges (Leong 2019). Another conventional laboratory method for estimating SWRC is to use the particle size distribution (PSD) curve and volume-mass properties (Wen et al. 2020). One of the theoretical methods is to apply the pore size distribution function (PoSD) and develop soil characteristic curve based on the soil pores network (Chen et al. 2019;Fredlund et al. 1994). It should be noted that, unlike PSD, PoSD is affected by climatic conditions (temperature, soil moisture content, number of rainy or rainless months) and even could change with soil shrinkage (Wen et al. 2021). Of course, some physical-based conceptual models have been propounded to predict SWRC. In these models, by using an appropriate PSD function and then conversion to PoSD, SWRC is modeled (Arya and Paris 1981;Assouline and Rouault 1997;Zhai et al. 2020). In most cases, particles and pores are assumed as solid spheres and void cylinders, respectively. However, the fundamental weakness raised about conceptual models is that they do not consider spheres in the aggregate or packed form in a specified volume. The present study focuses on this important gap so that by eliminating it an optimal function for PSD are provided. In recent years porous media structure prediction by computer modeling has become one of the researcher's favorite topics. Random close packing and pore network modeling (PNM) methods are used in many types of research and engineering sciences such as studying composite structures (Wu et al. 2010), distribution and loading of catalyst beds in chemical reactors (Dorai et al. 2012), study of different polymers (Ch et al. 2013), hydrocarbon-bearing reservoir rock (Al-Dhahli et al. 2013), microstructure and powder-based structures material (Benabbou et al. 2010;Zhou et al. 2009), packing cartons problems (Wu et al. 2017), etc. In mathematics, the problem of packing particles in a limited space is a subset of optimization problems. In the RCP method, solid particles are assumed as spherical, cylindrical, and cubic or pyramid (Li et al. 2010). Depending on the complexity and structure of considered media, using a combination of these figures is suggested. RCP method investigates the packing and relation of solid particles in scale under consideration. On the other hand, porosity, void spaces and channel connections are subject of the PNM method. In this method (for instance in 2D state), the channels and throats between particles are considered as circles and rectangles, respectively (Xiong et al. 2016). The Euler approach is the prominently used approach in RCP and the Lagrange approach is so in PNM method. Therefore, the results are not so different from each other. Particles could be used in the same size or with a size distribution (Baranau and Tallarek 2014;Daneyko et al. 2011;Santiso and Müller 2002). Packed media with spherical particles then compared on their porosity and coordination number. Particle density is defined as total solid volume to total volume ratio and coordination number is the number of particles in the vicinity of a specific particle (Anikeenko et al. 2008;Farr and Groot 2009). Porosity can be calculated from particle density. In the RCP method, mechanical contraction, Monte Carlo, particle drop and roll, spherical growth, and Lubachevsky-Stillinger (L-S) algorithm are among important rudimentary algorithms (basic algorithms) for the construction of porous media (Hitti and Bernacki 2013;Kansal et al. 2002;Kyrylyuk et al. 2009;Tobochnik and Chapin 1988). Hybrid or heuristic methods commonly are a combination of the above algorithms or algorithms with improved and coupled equations. For the mechanical contraction (sedimentation) algorithm introduced by Williams and Phillips (2003), first, all particles are generated in a finite space using a random distribution function. To minimize the space occupied by the particles and achieve the desired porosity, the particles' z-coordinate is reduced as much as possible based on their initial position. According to the definition of the non-overlap constraint for solid particles, each particle is sedimented vertically until it reaches the bed of stable stationary particles. The most important drawback of this method is the creation of unconventional spaces. This algorithm is strongly dependent on the initial position of the generated particles in the studied volume.
Although by changing the diameter of the particles and their initial position, this problem can be partially eliminated, but in any case, achieving the optimal arrangement is one of the main problems of this method. In the spherical growth algorithm, after producing the first particle as the central nucleus (core) of the set, considering the minimum particle distance and non-overlapping conditions, the first layer of particles is created in the vicinity of the desired particle. In this method, the growth of the layers surrounding the central core is circular (for 2D assumption). By increasing the circumfluent boundary of particles, the computation time increases significantly. Because the newly loaded particle sweeps the entire boundary of the set to achieve the optimal position (minimum distance from the central core) (Bargmann et al. 2018). He et al. (1999) proposed the Monte Carlo method for producing particles with different sizes in an arbitrary finite volume. According to the Monte Carlo algorithm, all particles are produced in completely initial random coordinates and smaller than their actual size. Then, by concerning the non-overlapping and volume interference constraints, the particles' volume is increased to achieve a dense arrangement with the desired density and porosity. Like the mechanical contraction method, this method is also affected by the initial position of the generated particles. However, in this method, by replication of the results (for example, up to 10,000 replication), the mentioned problem is eliminated and by averaging the results, the best arrangement and porosity are evaluated. For the drop and roll algorithm, developed by Visscher and Bolsterli (1972), unlike the Monte Carlo and the mechanical contraction methods in which all particles are generated simultaneously, the creation and evaluation of the stable position of the particles are investigated one by one. This method is more complex than previous methods and its calculation time is very long, But the results are more precise. It is worth noting that the RCP method is the starting point for research on porous media. Indeed, after modeling and generation of the structure of the porous media, various aspects such as the fluid flow, the rate of thermal resistance (or shape alteration), and the compressibility of the porous media are investigated. As mentioned before, the main objective of this study is to propose modified PSD and an upgraded mathematic model to predict coarse-grained soils structure as a typical granular porous media. The heuristic algorithm which is presented for the first time in this study is the combination of "spherical growth" and "dropping and rolling" by the implementation of modified functions in the matrix domain.

Modeling and system parameters
To predict the soil structure (sifted soil) RCP method is implemented by MATLAB coding. The heuristic algorithm applied in this study is the combination of "spherical growth" and "dropping and rolling" with some modified equations. Soil particles are assumed as spherical solid particles. In this model objective function is a random distribution of spherical particles in a finite cubic volume. Whereas soil media is continuum and continuum media modeling is difficult and more time-consuming, so in the first step, the mentioned cube (called finite volume) meshed up as depicted in Fig. 1. Discretizing and meshing precision is included by defining the resolution factor (RF). As shown in Fig. 1, just a set of specified nodes are utilized instead of the whole points of continuum media. These nodes represent their surrounding media. Obviously, more node leads to more precision which takes more run time and RAM capacity. Node numbers depend on the smallest particle size because the correct packing and occupation of each node by particles (spheres) must be guaranteed. Two functions to satisfy these conditions are defined as Eq. (1) and Eq. (2).
In the above equations, t cal points to calculation time, d i refers to particles (spheres) diameter, n i is the number of each particle, L c is cube edge size, and ε is the estimated porosity. In Fig. 1, grid spacing (called mesh-grid) is obtained by Eq. 3 as follows: In which L m notes two nodes distance and L c is cube edge size. Nodes coordination (x i , y i , z i ) is assumed as a center of spheres. It should be noted that in matrix space if the center of enough large sphere is allocated to an arbitrary grid (sphere whose radius is bigger than two neighbor node distance), necessarily some of the neighboring grids would be enclosed too. Therefore, necessarily each node is not at the center of a sphere. But surely every sphere is placed at the center in a randomly selected node (Fig. 2a). The next step is to define of a random function for the completely random placing of particles. In each computational step, a group of nodes occupied by each sphere with an optional radius which produced by the rand function in MATLAB. Normal distribution function, Weibull, and Lognormal are among the important random functions in 3-dimensional spaces. Many studies show the best approximation could be obtained using Weibull distribution for the size distribution of soils (Esmaeelnejad et al. 2016;Rhodes 2008). Therefore, Weibull and normal distribution functions are employed in MATLAB to construct spheres with random diameters. Neglecting these nodes for the next placement is obligatory since counterpart soil particles could not overlap each other (Fig. 2b). It's also an important point to reduce computational time.
As it's shown in Fig. 2b, in 2D space and for a general (i, j) node chosen as circle center point, neighboring nodes also placed inside of the circle so would not be considered again. This important condition is implemented by defining matrixes with zero and one arrays. The matrixes that include the number one array represent nodes that have not yet been occupied by the particle. For the next ball (with specified random radius) placement, the nearest node with minimum z to the last sphere would be selected according to Fig. 2a and c. Indeed, to increase the computational speed, the probability function is defined for the remaining nodes. To place the next sphere, the distance of all remaining nodes from the center of the previous sphere is calculated in the matrix domain. A node which is selected for the center of the next sphere must fulfill two constraints. The minimum distance from the center of the previous sphere and its minimum z i coordinate. The probability matrix is calculated for all remaining nodes and the node with the largest array (maximum probability) is selected. Obviously, in this case, there will be nodes with equal probability density (Fig. 2b). After assigning the node as the center of the sphere (with a random radius), the non-overlapping condition is tested by calculating the Euclidean distance of the node from the node allocated to the central core. If the above condition would not true, the node is stored in a matrix called NAN and removed from the calculation cycle, and the next node is chosen with the maximum probability. This cycle is continued until the particle loaded correctly and the next sphere loading calculations is resumed. As mentioned before, the present algorithm is a combination of "spherical growth" and "drop and rolling" algorithms by implementation in matrix space. To apply this idea, the cubic volume which contains particles is divided into m sublayers. Then by using conditional statements to constrain the model to necessarily fill the first layer with spheres and do not proceed to the next layer until the previous layer (or nodes) completely filled. After placement of each sphere, this sphere acts as a central (core) sphere for the next adjacent spheres. The 2-dimensional case of this process is represented in Fig. 2d. Sphere number represents its placing order. Flowchart of all calculations and particles placement in confined volume is shown in Table 1.
Our main research goal is to investigate contaminant dispersion in soil. In the first step, we need to have a good estimate of SWRC and SHCC. As mentioned before, one way to predict SWRC is to use the accurate PSD function. In the present study by applying a heuristic algorithm the structure and particle configuration of coarse-grained soils are investigated. Besides, properties such as porosity (or saturated volumetric humidity), saturated hydraulic conductivity and pressure drop are evaluated. These parameters are strongly influenced by soil texture and structure. Soil texture indicates the percentage of sand, silt, and clay particles and is an important criterion of soil PSD. However, the structure of the soil explains the configuration and arrangement of these particles. Therefore, we want to know this method of RCP to generate porous media which uses spheres, how closely resembles real media.

Materials and equipment
For preparing laboratory coarse-grained soils, a specified amount of soil is collected and sorted by different grain sizes using standard laboratory sieves which are used for gradation tests as it's shown in Table 2.
In order to ensure uniform PSD for each mesh, screening with a mechanical shaker is repeated three times. A plastic hammer was used during meshing to eliminate dust. Four sets with specific mesh sizes are shown in Fig. 3. Four types of coarse-grained soils are provided by mixing all of the above four particle size groups with mass weight mixing ratios as shown in Table 3. These mixing ratios are also used for modeling coarse-grained soils in MATLAB. For each size classification, Weibull and normal distribution functions are employed. Indeed, to generate spheres with random diameter, four Weibull and normal distribution functions are applied.
For a random variable X (diameter of spheres) that takes positive and continuous values, the Weibull probability density function is presented by Eq. (4) as follows: where k and λ refer to shape and scale parameters, respectively.

Porosity measurement
To porosity measurement of the mentioned laboratory coarsegrained soils, standard stainless steel cylinders (100.00 ml) were filled with them and saturated with water. Cylinders' weight was measured precisely before and after filling with soil and then water. Then the samples were placed in an oven at 105 °C for 24 h to ensure that all its moisture content evaporated. Samples' weight was measured again and porosity was calculated by its difference to previous measurements. For each soil, this experiment was repeated 3 times to ensure the accuracy of the results. The cylinders containing sifted soils are shown in Fig. 4.

Pressure drop and permeability
One important parameter that affects fluid flow in porous media, is pressure drop due to porous media resistance to flow through pores. Pressure drop increases by effective porosity reduction. A modified equation based on Kozeny-Carmen equation (Rhodes 2008) is utilized to investigate the pressure drop in the model. For a porous media with identical spherical grain, pressure drop could be expressed as Eq. (5): In which ΔP is pressure drop, L is length, µ is dynamic viscosity, v s is superficial velocity, ε is porosity, φ s is sphericity (which is equal one for a perfect sphere) and D p is equivalent grain diameter. Using Darcy law and rearranging Kozeny-Carman equation, the permeability of a porous media with identical sphere grain could be written as below: Permeability and saturated hydraulic conductivity relation could be expressed by Eq. (8): In which k is the permeability of porous media in m 2 , K s is hydraulic conductivity in m s −1 , µ is dynamic viscosity in kg m −1 s −1 , ρ is fluid density in kg m −3 and g is gravity acceleration in m s −2 (Loll et al. 1999). To precisely estimate the saturated hydraulic conductivity of groups of jammed sphere grains with different diameters, Irani and Callis (Irani and Callis 1963) equation is used (Eqs. 9a-9d). This applicable equation represents the equivalent mean diameter for a group of soil particles of different sizes. In which d g is the mean diameter of grains, σ g is the variance of grain size distribution, n is the number of different grain size groups, and f i is the mass percent of grains which has a diameter equal or less than M i in each group. Following the recommendation for application of geometric mean diameter to sieved sediments, M i is taken as the arithmetic mean of two consecutive particle size limits (Shirazi and Boersma 1984). Therefore, by using the mentioned equation and also using Weibull and normal distribution functions to generate spheres with different diameters, the mean diameter is calculated for four types of sifted soils.

Laboratory investigation of pressure drop
ASTM D 2434-68 is used for measuring pressure drop and saturated hydraulic conductivity in soil porous media. This standard is used to measure the pressure drop and hydraulic conductivity of disturbed granular soils in saturated conditions. The setup of the experiment is shown in Fig. 5. Sifted soils from previous steps (65 kg), are used for experiments.  Table 3 Mixing ratio for four soils by their mass weight a Type points out the assumed mixing ratios. For example, each 6 kg prepared soil type 1 consists of 3, 3, 1, and 1 kg of soil particles in the range of 0.3-0.6 mm, 0.6-1.18 mm, 1.18-2 mm, and 2-3.35 mm, respectively Type 1 a Type 2 Type 3 Type 4 0.3-0.6 mm 3 2 1 1 0.6-1.18 mm 3 2 1 1 1.18-2 mm 1 2 2 3 2-3.35 mm 1 1 2 3

Fig. 4 Prepared soil samples for porosity measurement
Environmental Earth Sciences (2021)

Experimental design method
By applying "Minitab" software and general full factorial design method, the order of experiments is completely randomized. Each experiment is repeated three times to assure the validity and accuracy of the results. The total number of 72 experiments are carried out. 6 treatments of 49, 59, 69, 79, 89, and 95 cm water head (constant water head) are used for four soil types. Due to using sifted soil particle which has specific grading, it is expected that pressure drop would be the same for the equal depth of column as soil prepared homogenously.

Modeling and laboratory results
The model result for type 1 soil structure is delineated as an example in Fig. 6. Also, the structure of the coarse-grained soil for the particle size distribution of 25− 3350 μm is shown in Fig. 7. Comparison of the results between normal and Weibull distribution functions showed that the Weibull distribution function has much better results in estimating soil porosity. In the next step, to reduce the difference between the estimated and empirical porosity, the skewness and elongation coefficients of the Weibull distribution function (in the terms of shape and scale parameters) were inclined to negative values. As a result, the use of the modified Weibull distribution function led to the development of an optimal RCP model for predicting porosity. According to the symmetric shape of the normal distribution function (around the mean value of each desired grain size range), its skewness coefficient is zero and its default elongation coefficient is 3. Comparison of the sifted soils porosity and the corresponding optimal RCP model showed that the average relative error between the two porosities is less than 4%. While for the normal distribution function, the predicted porosity showed larger values and the relative error percentage was up to 8.9%. The predicted values of porosity by mathematic optimal model and experimental data for four types of soils is presented in Fig. 8. The most important parameters in determining the porosity of the porous medium were replication, resolution factor, and particle size distribution function. Modeling of particles configuration carried out in two modes which were the least and most effective magnitude. While the rest of the parameters were assumed to be constant (see Table 4). The results of the sensitivity analysis showed that the porosity obtained from the model is strongly sensitive to the resolution factor and should be chosen with sufficiently large amounts. The results are represented in Fig. 9. By using a modified distribution function for sphere diameter and choosing a high-Resolution factor (RF > 250), the relative error of porosity was reduced to less than 3%. By normal distribution function and 50 < RF < 100, relative error up to 10% obtained. The resolution factor is a criterion of discretization accuracy and the number of nodes in the cube. For example, for a cube with a length of 1 cm and RF = 200, there would be 8 × 10 6 nodes. In this case, it's obvious that run time increases exponentially too (over 10 h for RF = 200 and a week for RF = 300). By approaching higher RF values and using the appropriate conditional clause in the Weibull distribution function, we could obtain spheres which closely approximate the real condition and build an optimum model with a more precise answer, but it needs higher computational power and run time. The predicted results for pressure drop reveal that Weibull distribution function gets a more accurate result than the normal distribution function (Fig. 10). As porosity value obtained by models overestimates real values, accordingly calculated pressure drops are less than measured   (Shirazi and Boersma 1984) proposed relation to calculate the arithmetic mean particle diameters for each grain classification. The average difference between actual pressure drop and estimated by the RCP model (using the normal distribution function) was 23.3%. While by applying the modified Weibull distribution function, this difference was reduced to less than 10%. As is expected, almost in all experiments such a result for ΔP/L was observed in manometers which is shown in Fig. 11. Saturated hydraulic conductivity of sifted soils (in the range of sandy soils) was measured using ASTM D 2432-68 and it was compared with model estimates. As an example, Estimated saturated hydraulic conductivity and measured values for soils type one and two are shown in Fig. 12. As previously explained, the modified Weibull distribution function is used to reduce the pressure drop difference to less than 10%. Subsequently, given Eqs. (5) and (7), saturated hydraulic conductivity values are predicted up to 10% larger than the actual values.

Conclusion
In this research innovative algorithm is employed to generate sandy soil porous media which is the combination of "dropping and rolling" and "spherical growth" algorithms.
The main advantage of this algorithm is the utilization of the continuous media discretization method and transforming to matrix space which could successfully model porous media. Also, by using applied equations in soil science, and a modified Weibull distribution function, an optimal RCP model was developed to predict the parameters of porosity, pressure drop, and saturated hydraulic conductivity. Comparison of the results of the optimal model and laboratory investigations reveals that the difference between the predicted and measured porosity is less than 4%. A difference of 10% is also observed for the parameters of pressure drop and saturated hydraulic conductivity. One of the unique characteristics of this model is building completely organized and structured porous media which happens for equal and monodispersed spheres in the range of 30 < RF < 70 (represented in Fig. 13).

Future study
Providing a suitable model to predict SWRC has always been a favorite of researchers. By applying a proper PSD and also ratios of fine sand/silty clay mixtures to their maximum densities, we can propose a conceptual model to prognosticate SWRC with acceptable accuracy. Indeed, after finding the PSD, using the suggested equations, we will try to obtain the corresponding PoSD function, and finally, SWRC will be predicted based on PoSD.