Understanding the characteristics of pore structures and their effects on permeability, elasticity, electrical properties, etc. is critical in a wide range of geoscience and petroleum engineering fields (Cai et al. 2014, 2020; Kerimov et al. 2018; Zhao et al. 2018; Wu et al. 2019; Xie et al. 2022; Zhu et al. 2022). Pore structure refers to the overall features of the size, shape, distribution, and connectivity of the pore space (Jiang et al. 2007). The size distribution and shape of particles as well as the binding modes between particles are important factors affecting the geometrical and topological properties of granular porous media (Cho et al. 2006; Baule et al. 2013; Hou et al. 2016; Kerimov et al. 2018; Li et al. 2022). These pore structures in turn affect the physical properties of porous media. Natural rock is one of the typical porous media, and it is challenging to theoretically investigate the effective physical properties of such materials due to the extremely complex and irregular pore structure.

Digital rock physics (DRP) is a flexible approach to study the physical properties of rocks (Andrä et al. 2013; Karimpouli and Tahmasebi 2016; Nie et al. 2019; Tan et al. 2021), which is based on imaging of pore space and component information, and numerically simulating of effective properties such as permeability, electrical resistivity, and elastic moduli. Dvorkin et al. (2011) and Andrä et al. (2013) summarized the DRP workflow into three steps: digital core modeling, pore space and mineral segmentation, and effective property simulation. Based on the DRP workflow, researchers have conducted extensive studies at the pore scale, including permeability (Karimpouli and Tahmasebi 2016; Miller et al. 2017; Yang et al. 2018; Liu et al. 2022), elasticity (Zheng et al. 2014; Karimpouli et al. 2018; Liang et al. 2020; Andhumoudine et al. 2021; Tan et al. 2021), and electrical properties (Liu et al. 2009a; Zhan et al. 2010; Zhao et al. 2013; Sun et al. 2014; Wu et al. 2020). Furthermore, with the rapid advances in computing and imaging techniques, the use of DRP procedures to study the physical properties is becoming a cost-effective avenue.

Digital porous medium modeling that is the first step in the DPR procedure is critical for determining the accuracy and breadth of petrophysical numerical simulations. In general, the methods for constructing digital models of pore and mineral structures can be grouped into two categories: physical experiment and numerical reconstruction (Blunt et al. 2013; Yang et al. 2023). The physical experiment approach is based on high-precision equipment such as optical microscope, nuclear magnetic resonance, and micro-CT to directly scan samples to gain digital descriptions (Madonna et al. 2013; Wu et al. 2019; Xia et al. 2019; You et al. 2021). However, the physical experiment method is time-consuming and costly, and does not have fully quantitative variability in modeling pore structures, thus affecting the breadth of DRP research. On the other hand, numerical reconstruction approach provides an inexpensive flexibility for building pore-scale images of porous media. For example, 3D pore images are stochastically reconstructed based on mathematical and statistical information from 2D scanned images, such as Gaussian field method (Adler et al. 1990), simulated annealing method (Talukdar and Torsaeter 2002), sequential indication method (Keehm et al. 2004), and multi-point statistics method (Okabe and Blunt 2004). Sandstone-like models are constructed by simulating the geological formation process of sedimentary rocks, the so-called process-based method (Bakke and Øren 1997; Øren and Bakke 2002). Random packing models of regularly or irregularly shaped particles are generated through the aggregation and transformation of particles, such as discrete element method (Pöschel and Schwager 2005; Al Ibrahim et al. 2019) and Markov process method (Tahmasebi and Sahimi 2018). Meantime, hybrid modeling approaches that combine several specific methods are becoming an increasingly important tool for generating pore structure models, such as CT and simulated annealing (Lin et al. 2017), process-based and simulated annealing (Liu et al. 2009b), Markov process and simulated annealing (Yao et al. 2013), and CT and mathematical morphology (Zha et al. 2017; Nie et al. 2019). Therefore, in this work, a hybrid method combining process-based and mathematical morphology is used to reconstruct 3D models of various pore structure types.

Changes in physical properties of porous materials usually follow regular trends, and obtaining these reference trends is of great significance for predicting physical properties and expanding applications (Andrä et al. 2013). In conventional laboratories, these trends are acquired by measuring lots of samples (Karimpouli and Tahmasebi 2016), which is not only time-consuming and costly, but also fails to elucidate the underlying mechanisms. However, DRP-based numerical simulations can be effectively used to compute physical property parameters, grasp variation trends, and explore the underlying mechanisms at the pore scale. Øren and Bakke (2003) explored the permeability, conductivity, capillary pressure and other properties of 3D porous medium models based on the numerical simulation method. Wojtacki et al. (2017) and Nie et al. (2019) revealed the changing laws of elastic moduli and permeability by uniformly transforming the pore space of sandstone and carbonate rocks, respectively. Han et al. (2020) discussed the relationship between pore structure and elastic properties of granular rocks by defining a new parameter cementation exponent. Zhao et al. (2021) explored the acoustic feature and feasible classification criteria of carbonate reservoirs with different pore structure types. Although the DRP workflow enables numerical studies at the pore scale, the relationship between pore structure and effective properties is still not understood and described accurately (Han et al. 2020). This can be attributed to two factors. On the one hand, the response mechanisms and sensitivities of different rock physical properties to the pore structure are distinct dramatically, that is, there are significant differences in the variation laws of pore structure on different physical properties. On the other hand, the pore structure of porous media has a wide range of complexity, which makes the influence trends often difficult to be described using specific models.

Due to the aforementioned issues, in this paper, a workflow is improved to specifically study the effect of pore structure on physical properties. It focuses on exploring the micro-scale mechanism of the pore structure-physical properties relations, and an example is performed based on the numerical reconstruction models. Firstly, the hybrid approach combining process-based and morphology-based methods is described, and several groups of pore structure models are constructed. These models serve as the basis for discussing the relationship between pore structure and physical properties and analyzing their mechanisms. Then, the pore structure features and differences between these models are compared, which is an important step in understanding the pore structure-petrophysics relations. Finally, the variation laws of the influence of pore structure on physical properties (absolute permeability and elastic modulus) are discussed, and the mechanisms are explained at pore scale.

With reference to standard DRP concept where rock physical property trends are driven by a large number of samples (Karimpouli and Tahmasebi 2016), a workflow for studying effective physical properties based on extensive pore structure models is presented. Figure 1 illustrates the three steps of this workflow: (1) generating the digital models with multiple types of pore structures, (2) analyzing the characteristics and differences of the model pore structures, and (3) computing the volume effective physical parameters, and discussing the relationship between the pore structure and the effective properties and the microscale mechanism.

Workflow for studying effective physical properties based on extensive pore structure models. a Constructing lots of 3D pore structure models; b Analyzing pore structure characteristics; c Calculating physical parameters, and analyzing trends and mechanisms

Compared to standard DRP procedure, this workflow focuses more on pore structure characteristics and their impact on effective properties. In the first step (digital modeling), the focus is on considering the diversity of pore structures. In the second step, the point is not on achieving precise segmentation of multi-components, but on analyzing the characteristics and differences of pore structures. The focus of the third step is not to discover the various factors affecting the physical properties, but to analyze the laws and mechanisms of the effect of single factor (pore structure) on the effective parameters.

2.1 Pore structure modeling

To illustrate the workflow of exploring the relationship between pore structure and petrophysical properties and microscopic mechanisms, it is first necessary to establish the models with specific evolutionary characteristics on the pore structure. A hybrid approach combining process-based and morphology-based techniques to reconstruct pore structure models is presented. The grain packing formed by the process-based is used as the initial model, and morphological operations are then performed to transform the pore space to generate models of various types of pore structures.

2.1.1 Process-based method

In the complex geological environment and long-term diagenesis, the formation and evolution of rocks usually undergo extremely complex physical, chemical, and biological processes. The researchers propose a method for reconstructing porous rocks by simulating the formation process of sedimentary rocks, namely process-based approach, which can be summarized into three basic procedures: sedimentation, compaction, and diagenesis (Bakke and Øren 1997; Øren and Bakke 2002).

The particle size distribution affects the pore structure properties of the model results during sedimentation (Li et al. 2019; Luo et al. 2020). The process is simulated as settling and accumulation of clastic particles, reaching a steady state. Grain packing that has a large porosity is the basis for the subsequent transformation of pore space by mathematical morphology. In addition, digital core models similar to sedimentary rocks can be formed through the simulation of compaction and diagenetic processes. The entire modeling process can be quantitatively controlled, and the reconstructed models resemble the actual cores in many features such as porosity, geometry, and connectivity (Biswal et al. 1999; Øren and Bakke 2002).

2.1.2 Mathematical morphology

Mathematical morphology has been widely used in 2D or 3D image processing. The method consists of a series of operations, including the four basic algorithms of dilation, erosion, opening and closing, as well as other advanced algorithms that are combinations of basic algorithms (Serra 1982; Hilpert and Miller 2001; Liu et al. 2009a). These algorithms can be used to implement a variety of image processing, including image classification and identification, raw image reconstruction, feature extraction, pore-throat identification (Silin et al. 2003; Xu et al. 2020), and pore space reconstruction (Zha et al. 2017; Nie et al. 2019).

In the field of digital rock physics, fluid distribution is usually simulated by the open operation (Hilpert and Miller 2001; Liu et al. 2009a; Sun et al. 2014; Andhumoudine et al. 2021; Tan et al. 2021). Referring to similar means, the pore space can be reconstructed by applying a series of morphological operations to generate diverse pore structure models. The improvement of this modeling approach over fluid distribution simulation mainly includes the application of various morphological algorithms, not just the opening, and in addition, the assumption that matrix cement grows in the pore space to form specific pore structures.

The erosion and dilation algorithms (Liu et al. 2009a) are defined as

where A is the initial set, B is the structuring element, and (B)_x represents the translation of B by x.

The opening algorithm is represented as the set A is eroded by the structuring element B, and then the result is dilated by B, while the closing algorithm does the opposite. The opening and closing algorithms are represented as

To illustrate the morphology-based pore structure processing, the basic algorithms are used to process 2D porous media images as an example (Fig. 2).In the figure, black indicates the solid matrix, white indicates the pore space, and gray represents the processed part. Figure 2a is the original image, and the results obtained by dilation, erosion, and opening algorithms are shown in Figs. 2b–d, respectively. The distribution and percentage of solid matrix (or pore space) can be altered by morphological operations, which provides feasibility for morphological-based modeling.

Schematic of 2D porous image processing using basic morphological algorithms. a Initial image; b Dilation result; c Erosion result; d Opening result

2.2 Calculation of rock physical properties

The calculation of rock physical properties is a considerable part in the DRP workflow. In this work, two important properties, permeability and elasticity (including elastic modulus and wave velocity) are investigated. Various algorithms for digital rock property simulation are summarized by Andrä et al. (2013).

2.2.1 Quasi-static pore network model

The methods of permeability calculation of 3D porous images are usually summarized into two categories. One is to directly simulate the flow field based on digital images using lattice Boltzmann or finite difference methods (Miller et al. 2017; Yang et al. 2018; Liu et al. 2022), and then calculate the permeability parameters. The other is to establish the pore network by extracting the pore structure information from digital image, and then calculate the flow parameters of the model based on quasi-static conditions (Dong and Blunt 2009; Alyafei et al. 2016).

The quasi-static pore network model is an efficacious approach for simulating flow processes and estimating permeability parameters. The first step is to extract the pore network that is equivalent to the digital core geometrically and topologically. Then the pressure difference is applied between the inlet and outlet of the network model based on the quasi-static condition, and the capillary pressure and the flow rate of each phase fluid are measured. Finally, the permeability parameters are calculated by semi-analytical formula (Øren and Bakke 2002; Ryazanov et al. 2009; Alyafei et al. 2016). The semi-analytical approach is flexible for multiphase flow determinations where capillary forces dominates, so it saves time for calculations on large numbers of digital cores. The extraction of the pore network facilitates the analysis of the pore structure features of the digital core images, which is a powerful tool for further exploration of the relationship between structural parameters and rock physical properties.

2.2.2 Static elastic finite-element method

The finite element method proposed by Garboczi and Day (1995) is used to numerically simulate the elastic properties of digital porous media. The model needs to be divided into non-overlapping grid units during the calculation. The digitized binary image satisfies this condition well, and each pixel is marked as an element unit. Macroscopic strains are applied in different directions in which the material is stressed, and the elastic potential energy is calculated as a quadratic polynomial of the elastic displacement on each pixel unit. The stress and strain values are calculated by solving the minimum elastic strain energy using the conjugate gradient algorithm (Arns et al. 2002; Jouini et al. 2015; Andhumoudine et al. 2021). Based on the calculated stress and strain components (τ_ij, ε_ij, where i, j = 1, 2, 3), the effective bulk modulus (K) and shear modulus (µ) can be obtained by the following equations (Garboczi and Day 1995):

Using the calculated elastic moduli and the known material density, the P- and S-wave velocities are determined by

A hybrid approach combining process-based with mathematical morphology is proposed to generate models with various pore structure types, and the graphical flowchart is illustrated in Fig. 3. First, process-based method is applied to produce grain packing models with diverse size distributions and shapes that have stable contact states through particle settling one by one (Fig. 3a). Then, the grain packing is binarized to provide an initial model for morphology-based pore space processing (Fig. 3b). The initial model is processed using the morphological algorithms so that its pore space can be directionally reconstructed to form new pore structure models, as shown in Fig. 3c.

Hybrid modeling workflow combining the process-based and morphology-based methods. a Sedimentary process simulation; b 3D digital image processed by binarization; and c Pore space reconstruction based on mathematical morphology

The mathematical morphology method is introduced into the digital core modeling by simulating the cementation process of sedimentary rocks. The purpose is to achieve the evolution of multiple cementation modes and establish models with gradual changes in pore structures. Three representative pore structure models simulating sandstone cementation are constructed. The first is to apply the dilation algorithm to realize uniform expansion of the surface of solid particles, which is similar to the simulation of uniform cementation in process-based modeling. The second is to preferentially treat the smaller pore space using the closing operation, allowing the pore space to be converted into solid matrix. The third is to use a combination of algorithms (ie, assuming that the entire 3D image is 1, then subtracting the result of closing operation, and adding the initial model) to preferentially transform larger pore space. Figure 3c shows the results of the directional transformation of the pore space using different morphological algorithms, where black represents the solid matrix, white represents the pore space, and red represents the new solid matrix formed by reconstruction. It is worth noting that in addition to the three algorithms applied here, other operations or combinations of operations can be performed to transform binary models according to specific needs.

In pore network representation, pore space is represented as pore network consisting of pores and throats, where pores correspond to larger void spaces, whereas throats correspond to narrow regions connecting pores (Dong and Blunt 2009). In order to facilitate the comparative analysis of models, the above three types of models are named as uniform reconstruction model, throat reconstruction model, and pore reconstruction model based on the concept of pore network and the characteristics of morphological reconstruction. These models realize the directional transformation in the pore space structure, thus greatly altering the porosity, geometry, and topological properties of the initial model, which ultimately affect the petrophysical properties. As for the calculation time of the reconstructed model based on morphology, it is comprehensively affected by morphology algorithm, structure element size, and image size. Taking the image size of 400³ pixel as an example, it takes 8 min to perform the opening operation with the structure element size x = 5, and several hours to establish the 15 submodels used in the work.

This hybrid modeling approach combines the advantages of process-based and mathematical morphology methods, including: (1) satisfying the stable packing of grains under gravity constraints, and the initial model has a large porosity; (2) quantitatively adjusting the size distribution and shape of grains; (3) directionally transforming the pore space by applying diverse morphological algorithms; (4) progressively changing the porosity, geometric and topological properties; (5) realizing faster modeling speed and time-saving.

As matched with an improved workflow, this hybrid modeling approach focuses more on transforming pore spaces to generate models of multiple pore structure types, corresponding to the first step in the improved workflow. In addition, traditional workflows focus on threshold segmentation in this process. Unlike physical modeling method, this method does not need to focus on threshold segmentation of different components in the workflow, as different components in numerical modeling method are predetermined. Furthermore, numerical modeling differs from physical modeling in that it is limited by instrument scanning resolution and can accurately define skeletons and pores. Therefore, the focus of the second step in the improved workflow is not on achieving precise segmentation of pore space, but on analyzing the characteristics and differences of pore structure.

4.1 Pore structure models

4.1.1 Grain packing model

Changing the particle size distribution affects the porosity and pore structure of the model to some extent (Lee and Lee 2013; Kerimov et al. 2018). The models of varying sized particles can be obtained by process-based simulation. In terms of particle size characterization of rocks, sorting is an effective index for the overall evaluation, and the sorting coefficient S is one of the commonly used parameters for characterizing particle distribution properties. The sorting coefficient is expressed as the ratio of the standard deviation to the average value for all particle (Sohn and Moreland 1968). The smaller S is, the more concentrated the grain distribution is, and vice versa, which is defined as

where, σ is the standard deviation of grain size distribution; E is the average.

A set of grain packing models are constructed by varying grain size distributions. To make it easier to observe changes in the model pore structures, uniform distribution functions are selected for the deposition grains, so that the grain size distribution becomes wider as the range of the distribution function expands. This makes the porosity and pore structure change progressively in the model results.

4.1.2 Morphology-based reconstruction model

Porous media often produce different pore structure in the process of formation or establishment. For example, various kinds of cementation and diagenesis occur that change the pore structure with the change of the geological environment during rock formation. The modeling workflow is detailed in the previous section, and the following are some specific examples of applying the method.

Mathematical morphological operations are performed on the digital packing model to directionally transform the pore space, thereby generating new pore structure models (Fig. 3c). The sub-model No. 10 in the packing model is used as the initial model for morphological modeling to ensure the consistency of the initial pore structure. Models of different pore structure types are constructed by applying the three morphological algorithms described above. In addition, changing the radius of the structuring elements enables the porosity to vary gradually. Thus, a significant number of 3D digital core models with gradually changing pore structures are finally obtained.

4.2 Pore structure characteristics

The geometry and topological structure of porous media affect the physical properties. One of the important purposes of establishing rock structure models is to discover the possible trends of petrophysical properties. This hybrid modeling approach flexibly build multiple sets of models with significant differences in pore structure. It is easier to explore the mechanisms and regularities of rock physical properties at the pore scale by analyzing the pore structure characteristics of these models.

Four types of models, including grain packing and morphological-based reconstruction models, all contain 15 sub-models, and their porosity evolution is shown in Fig. 4. Overall, the porosity of each type of model is gradually changing. There is a slight decrease in the porosity of the grain packing model, however, the porosity of the three kind of models of morphology-based shows a significant decrease because the pore space is transformed greatly when the morphological algorithms are applied. The porosity and pore structure vary over a wide range, providing a model basis for investigating the physical properties of porous media.

Porosity evolution of four kinds of pore structure models. The porosity of the grain packing model has a small variation range, while the morphological-based reconstruction models have a large porosity variation

4.2.1 Grain packing model

Figure 5 shows the 3D representation of the three sub-models in the grain packing model, including their packing images, binary images, and pore networks. Generally, the pore space of the sub-model changes continuously as the grain size distribution function changes. The sub-model (No. 1) only contains grains of one size, so the grains are arranged uniformly, and the pore structure are relatively regular (the first column in Fig. 5). The sub-model (No. 8) contains grains with a wider size distribution, resulting in tighter grain packing and smaller pore spaces. As the distribution function expands the sub-model contain more small particles, and therefore generally have tighter packing and more complex pore structures (sub-model No. 15 in Fig. 5). The changing characteristics of these pore network models correspond to the pore structure analysis, as shown in Fig. 5c. Meanwhile, the porosity of the grain packing model has a slightly decreasing trend (Fig. 4), which is consistent with the grain model results (Kerimov et al. 2018; Al Ibrahim et al. 2019).

3D representation of three sub-models in the grain packing model, including a Packing image, b Binary image, and c Pore network

4.2.2 Morphology-based reconstruction model

The morphology-based reconstruction model greatly changes the characteristics of the original model such as porosity, grain contact form, and pore connectivity. The characteristics and differences of the three kinds of models can be grasped by comparing their pore structures. When observing and comparing three types of models, it is necessary to consider the evolution process simulated during rock modeling. Therefore, it is important to observe the changes of these models from the direction of decreasing porosity in order to facilitate the understanding of pore structure evolution.

For each type of model, five sub-models (with porosity close to 30%, 25%, 20%, 15% and 10%) are selected for visual comparison of binary image (Fig. 6) and pore network (Fig. 7), which is conducive to the comprehensive understanding of the variation characteristics in the pore structures. For the uniform reconstruction model (Figs. 6a and 7a), as the particle surface grows uniformly, the pore space is continuously compressed, thus changing the porosity and pore structure of the initial model. As the porosity decreases, the number and size of pores and throats, as well as the coordination number, gradually decrease (Fig. 7a), resulting in a continuous weakening of connectivity.

Comparison of binary 3D images for three types of pore structure models. a Uniform reconstruction model; b Throat reconstruction model; c Pore reconstruction model

Comparison of pore network models for three types of pore structure models. a Uniform reconstruction model; b Throat reconstruction model; c Pore reconstruction model

For the throat reconstruction model, the smallest pore space in this model is preferentially transformed into solid matrix, resulting in a more concentrated pore space (Fig. 6b). The number of throats in the pore network models decreases dramatically, so that the connectivity between pores suddenly weakens (Fig. 7b). When the porosity approaches 10%, the pores still remain abundant, but the throats are reduced to almost zero and the pores are completely isolated, which results in the complete obstruction on fluid flow. On the other hand, this form of pore structure increases the contact area between particles, which is beneficial for enhancing the stiffness of the solid matrix, and ultimately making the model more resistant to compression.

It is also essential to compare pore reconstruction model that preferentially remodel along the large pore space, which is filled into solid matrix with successive morphology-based processing (Fig. 6c). The pore space is continuously compressed and becomes more dispersed. Therefore, the number of pores and throats increases in the intermediate porosity. In terms of elasticity, the contact points between solid particles in this model do not form strong cementitious matrix, resulting in weaker compressive strength than other models.

The geometric parameters of three types of models, including total porosity, effective porosity, and specific surface area, are calculated (Fig. 8). Figure 8a shows the variation of effective porosity versus total porosity. The solid diagonal line indicates that the effective porosity is equal to the total porosity, and the deviation between the data points and the diagonal line reflects the change in the proportion of effective porosity. The feature of the throat reconstruction model is that the narrow pore space of the original model is reformed and cemented preferentially, so when the total porosity is large (about 0.28), the effective porosity deviates from the diagonal line and decreases rapidly. When the porosity is about 0.18, there are no effectively connected pores at all. The effective porosity of the uniform and pore reconstruction models is basically consistent with the diagonal, indicating that the pore spaces of these models have good connectivity. Their effective porosity does not drop to zero until the total porosity is around 0.05, and this porosity feature is consistent with the limiting connected porosity (Schwartz and Kimminau 1987).

Comparison of three types of pore structure models in effective porosity and specific surface area versus total porosity. a Effective porosity; b Specific surface area

Furthermore, Fig. 8b shows the variation in pore specific surface area versus total porosity. Specific surface area is described as the ratio of total surface area to total volume in pore space, and its unit is μm²/μm³. As the porosity decreases, the pore space of the three types of models is quite different, so the specific surface area results are significantly distinct. Since the large pore space is preferentially processed in the pore reconstruction model, the large pores are partitioned into more small pores, resulting in a larger value of the specific surface area than the other two types of models. However, the throat reconstruction model has the smallest variation, and the uniform reconstruction model is between them. The comparison of effective porosity and specific surface area illustrates that there are significant differences in the pore structures among the three models.

The pore networks of the models are extracted using the maximum ball algorithm (Dong and Blunt 2009), and the pore throat parameters are counted. Figure 9 shows the comparison of network parameters extracted from three types of models, including the number of pores and throats, the average radius of pores and throats, and the average coordination number.

Comparison of three types of pore structure models in network parameters. The parameters include a the number of the pores, b the number of the throats, c the average radius of the pores, d the average radius of the throat, and e average coordination number

The variation curves of the total number of pores and throats are shown in Figs. 9a, b, respectively. In general, there are great differences among the curves. The pore reconstruction model is basically the largest, followed by the uniform reconstruction model, and the throat reconstruction model is the smallest in terms of the number of pores and throats when the porosity is the same. The throat number of the throat reconstruction model decreases rapidly when the porosity is still large (Fig. 9b), and this change profoundly reflects the modeling characteristics of the closing operation preferentially cementing the narrow pore space (throat).

The variation curves of the average radius of pores and throats are indicated in Figs. 9c, d, respectively. The mean radii of the uniform and pore reconstruction models gradually decreases with decreasing porosity. However, the average radii of the throat reconstruction model have an upward trend, which is because the small pore space is preferentially transformed. It should be noted that although the average radii of pores and throats in the throat reconstruction model is the largest, it does not imply that the model has higher permeability. On the contrary, its permeability decreases much faster than other models (the effect of pore structure on permeability is discussed in the next section).

The coordination number reflects the connectivity between the pores and is one of the key factors affecting permeability. Figure 9e illustrates the change in the mean coordination number for the three models. As the porosity decreases, the coordination number of all models decreases gradually, but the decreasing trend is different. When the porosity is greater than 0.2, the uniform and pore reconstruction models decrease slowly, while the throat reconstruction model drops steeply.

Overall, there are significant differences in the evolution of pore structure among these three models, and it is these characteristic differences that are conducive to exploring the mechanisms and laws of the influence of pore structure on physical properties of materials. In traditional workflows, the emphasis of the third step is to analyze various factors affecting petrophysical properties, while the improved process focuses on discussing the relationship between individual factors (pore structure) and effective parameters, as well as the mechanism analysis.

4.3 Permeability of pore structure models

4.3.1 Grain packing model

The particle size distribution affects the size, shape, connectivity, etc. of the pore space, and thus has a remarkable impact on permeability (Lee and Lee 2013; Kerimov et al. 2018). The relationship between pore structure parameters and permeability is investigated based on the grain packing model, as shown in Fig. 10. The permeability of the model increases with the increase of the average grain radius, and this relationship can be well fitted using exponential formula (Fig. 10a). The relationship between porosity and permeability can be fitted by exponential equation, as shown in Fig. 10b, which provides a reference for the changing trend of permeability of porous media with larger porosity. In addition, the relations of sorting coefficient and specific surface area versus permeability can also be fitted using exponential equations, as shown in Figs. 10c, d, respectively. The coefficients of determination for all fitted curves are greater than 0.97, indicating that there are fairly close relations between permeability and pore structure parameters in the grain packing model. The variation trend of permeability is consistent with the pore structure analysis. Particle-filled porous materials are widely used in geosciences, chemical engineering, and other fields, so the exploration of these relationships provides a theoretical reference for understanding permeability laws.

The relationships of various pore structure parameters versus permeability in the grain packing model. The pore parameters include a Mean grain radius, b Porosity, c Sorting coefficient, and d Specific surface area

4.3.2 Morphology-based reconstruction model

The numerical results of permeability for the three types of morphology-based reconstruction models are shown in Fig. 11. Overall, the permeability of all models declines with the reduction of porosity, but the trends are markedly different. This is because there are considerable differences in pore structure evolution of the models, which have a significant impact on permeability.

Variation trends of porosity and specific surface area versus permeability for three kinds of pore structure models. a Uniform reconstruction model; b Throat reconstruction model; c Pore reconstruction model

The relationship between porosity and permeability in the uniform reconstruction model can be fitted by power formula (Fig. 11a), which illustrates that the permeability has a considerable downward trend as the porosity decreases. When the porosity drops to 5%, the permeability decreases to zero, indicating that this is the limit state of fluid flow in the pore space, and the limit porosity is consistent with the results (Schwartz and Kimminau 1987). In addition, Fig. 11a also shows the exponential changes in specific surface area versus permeability. The permeability gradually decreases with the reduction of specific surface area, indicating that the uniform compression of pore space has a significant impact on the permeability.

The permeability of the throat reconstruction model is calculated (Fig. 11b). When the porosity is greater than 0.35, the decreasing trend of permeability is smaller than that of uniform and pore reconstruction models. There is a steep decline in the range of porosity from 0.20 to 0.35, and the value drops to zero although at really high porosity (0.20). This changing trend of permeability is controlled by the variation characteristics of pore structure. The smallest pore spaces in this type of model are transformed initially, which mainly are small pores or throats that have a slight contribution to permeability, so permeability decreases slowly in the initial stages. As the larger pore spaces are processed as solid matrix, the connectivity of the model is greatly weakened because these pore spaces are channels for fluid flow, causing a sharp drop in permeability at this stage. The parameters of porosity and specific surface area have no clear relations with permeability in this type of model.

The permeability changes of the pore reconstruction model are shown in Fig. 11c. The permeability decreases gradually with the reduction of porosity, and their relationship can be well fitted by applying power formula. The limit porosity of pore space seepage is also close to 5%. Simultaneously, the exponential equation can be applied to perfectly fit the relation of specific surface area versus permeability (the coefficient of determination is 0.99). The variation trend and fitting relation of this type of model are tremendously different from those of the throat reconstruction model, which is attributed to the obvious discrepancy in the variation characteristics of pore space.

Furthermore, the permeability results of the three models are displayed together for comparison, which is beneficial to observing the changes in the effect of pore structure differences on permeability. Figure 12 shows the comparison of changes in porosity versus permeability among the three models. The ordinate is set to linear and exponential forms. The detailed analysis of each type of model has been explained above, and here is just a visual comparison. In general, the permeability of the three models varies gradually with the decrease of porosity, but the change patterns are diverse. When the porosity is in the range from the maximum to 0.3, the permeability of the throat reconstruction model is higher than that of the other two models. This is because the narrow pore spaces in this model are transformed first, yet their contribution to fluid transport is relatively low, which is agreeable with the understanding that flow channel size has a quadratic response to permeability in the capillary model. In the porosity interval of 0.2–0.3, the permeability of the throat reconstruction model drops steeply with the reduction of porosity, whereas there is a power function decreasing trend for the other two kinds of models (Fig. 11). This is attributed to the fact that the transport channels of the throat reconstruction model are rapidly transformed into solid matrix, while the other two types of models have extensive flow channels (Fig. 7). When the porosity is less than 0.2, the downward trend in permeability for the uniform and pore reconstruction models becomes more pronounced with decreasing porosity (Fig. 12b).

Comparison of changes in permeability for three kinds of pore structure models. The vertical coordinate values vary a Linearly and b Exponentially

It can be seen from the comparison that the three models with different pore structure evolution have distinct variation trends in permeability, and the discrepancy in the pore structure is the main factor causing permeability results. The possibility of permeability variation is explored based on these comparative analysis, and some relations between pore structure parameters and permeability that cannot be quantitatively measured in the laboratory are established (Fig. 11).

4.4 Elasticity of pore structure models

Unlike permeability analysis, which primarily focuses on pore space connectivity (topology), elastic properties analysis takes more consideration into solid phase characteristics such as solid composition, cement content and type, particle contact mode, etc. Since the main purpose of this work is to discuss the influence mechanism and trend of pore structure on elasticity, it is assumed that all pore structure models should be set with the same mineral composition during the numerical simulation. Specifically, it is assumed that the solid matrix and the pore space are quartz and water, respectively, and the physical parameters of the mineral components (Wei et al. 2018) are shown in Table 1. Therefore, the relationship between pore structure and elastic properties and the micro-scale mechanism can be studied.

4.4.1 Grain packing model

The grain packing model is built by varying the size distribution, which affects the compactness of particle filling, so it has an impact on the elasticity (Madadi and Saadatfar 2017; Kerimov et al. 2018). The effective elastic bulk and shear moduli are calculated numerically using the finite element method, and the relationships between pore parameters and elastic modulus are plotted (Fig. 13). These relationships can be fitted using linear formulas, which are shown in the Fig. 13. Both bulk modulus and shear modulus decrease steadily with the increase of mean grain radius and porosity, but increase with the increase of sorting coefficient and specific surface area. Moreover, the coefficient of determination (R²) are more than 0.91, indicating that the pore structure parameters are significantly related to the elastic properties in the packing model.

The relationships of various pore structure parameters versus elastic modulus in the grain packing model. The pore parameters include a Mean grain radius, b Porosity, c Sorting coefficient, and d Specific surface area

4.4.2 Morphology-based reconstruction model

In practice, even though porous media have the identical porosity and mineral composition, their elastic modulus (or wave velocity) properties are sometimes different, which is mainly due to the differences in the grain contact patterns (Sun et al. 2021). For example, various complex actions control the size, shape, and cementation patterns of solid particles during rock formation, which greatly affects the elastic characteristics of rocks. Therefore, it is important to explore how and to what extent the pore structures affect the elastic modulus or wave velocity based on digital rock physics technology.

Figure 14 shows the comparison of the effective elastic modulus and wave velocity results for three types of morphology-based reconstruction models. Overall, the elastic moduli and wave velocities of all models rise significantly with decreasing porosity, again proving that porosity is one of the most important factors affecting elastic properties. On the other hand, the variation curves of these models have different evolutionary laws, which is caused by the huge difference in pore structure.

Comparison of changes in elastic parameters (elastic moduli and wave velocities) for three kinds of pore structure models. a Bulk modulus; b Shear modulus; c P-wave velocity; d S- wave velocity

The blue dots in Fig. 14 represent the elastic results of the throat reconstruction model. The elastic moduli and wave velocities of this model are higher than those of the other two models when their porosity is the same. This can be explained that the pore structure of the model is characterized by preferentially transforming the narrow throats of the initial model, so that the grain contact points are cemented, and the solid matrix tends to be combined as a whole (Fig. 6b). As a result, the material formed by this cementitious structure has stronger compression resistance and greater stiffness. The elastic parameters of this model are increased obviously and are higher than those of the other two types of models in the initial stage (i.e., the porosity is close to 0.35). This can be explained by the fact that the overall stiffness of the material can be significantly improved as long as there are very few parts of matrix cemented at the contact points between the grains. In contrast, the elastic parameters of the pore reconstruction model indicated by the green dots are generally smaller than those of the other two models (Fig. 14). This is attributed to the changing characteristics of the pore structure of this model, where the large pore spaces (pores) are remodeled first to make the growth of cementitious matrix away from the particle contact points (Fig. 6c). Materials with this pore structure feature has relatively weak compression resistance and low stiffness. In addition, the variation curve of the uniform reconstruction model is located in the middle (red dots in Fig. 14). The relations of elastic moduli and wave velocities versus porosity varies approximately linearly, while the other two models are nonlinear. This can be understood that the pore structure of this model is uniformly changed, which only leads to a steady decrease in porosity, but no significant changes in pore structure characteristics (or grain contact patterns) (Fig. 6a).

The variation trends of the elastic parameters of the three types of models are further compared. When the porosity is maximum, the elastic moduli and wave velocities of the three models are the same due to the same initial model used. As the porosity falls slightly, the manifestation of the effect of pore structure differences on elastic properties appears rapidly. When the porosity decreases to 0.3, the difference of elastic moduli between the throat reconstruction model and the pore reconstruction model reaches the maximum, the wave velocity difference at this time is also the largest, indicating that there is an evolution period that makes the influence of the pore structure on the elasticity to the maximum. Specifically, comparing the difference between throat reconstruction and pore reconstruction models, the bulk modulus is 4.08 GPa (13.81–17.89), and shear modulus is 5.31 GPa (10.33–15.64). The P-wave velocity is 0.68 km/s (3.59–4.27) and the S-wave velocity is 0.51 km/s (2.20–2.71). Furthermore, the change rates of these four elastic parameters reached 29.54% (4.08/13.81), 51.40% (5.31/10.33), 18.94% (0.68/3.59), and 23.18% (0.51/2.20), respectively. From the data comparison, it can be further found that the pore structure has a greater influence on the shear modulus than the volume modulus, and the S-wave velocity than the P-wave velocity. The differences in the elastic results of the three models decrease gradually with the further reduction of porosity, that is, the three curves converge gradually, indicating that the effect of pore structure difference on the elasticity has a weakening trend. When the porosity is less than 0.1, the contrast of elastic parameters between these models becomes small. This is because even though the difference in pore structure is still large at this time, the number of voxels where the pore structures of the three types of models differ have been dramatically reduced due to the small porosity.

The evolution law and micro mechanism of the effect of porosity and pore structure on the elastic properties are discussed by numerical simulation. In conclusion, porosity is one of the most important factors affecting the elastic properties, and the pore structure also greatly affects the elastic properties. The pore structure has a certain degree of influence on the elasticity, but the degree is weaker than that of the porosity. Numerical simulations based on this workflow provide insights into the mechanisms and trends in the effects of pore structure features on the elastic properties of porous media. These studies reveal that pore structure evolution is extremely complex, but the relationship between macroscopic properties and microscopic pore structures can be better understood through the improved workflow analysis.

This paper focuses on a specific workflow for studying the effects of univariate pore structure on the petrophysical properties. Several types of pore structure models are reconstructed based on the hybrid modeling method, and the relations between pore structure and permeability/elasticity and the microscopic mechanisms are discussed.

1. (1)

   Improving a workflow that focuses on discussing the relations between pore structure and petrophysical properties, as well as pore scale mechanisms.

2. (2)

   A hybrid method combining process-based and morphology-based is realized, which can transform pore space specifically to simulate the evolution of rock cementation.

3. (3)

   The relationship models between pore structure parameters and permeability/elastic parameters in the grain packing model are established. The effect of pore structure evolution on permeability/elastic parameters and the microscopic mechanism in three types of morphology-based reconstruction models are analyzed.

4. (4)

   The influence degree of pore structure on elastic parameters (bulk modulus, shear modulus, P-wave velocity, and S-wave velocity) is quantified, reaching 29.54%, 51.40%, 18.94%, and 23.18%, respectively.