Estimating Permeability of Porous Media from 2D Digital Images

: Digital rock physics (DRP) has been widely used as an effective approach for estimating the permeability of porous media. However, conventional implementation of DRP requires the reconstruction of three-dimensional (3D) pore networks, which suffer from intensive memory and underlying uncertainties. Therefore, it is highly signiﬁcant to develop an approach only based on two-dimensional (2D) cross-sections of parent samples without 3D reconstruction. In this study, we present a novel approach that combines the Kozeny–Carman equation with fractal theory to derive a bridge function that links 2D cross-sectional images and 3D pore structures of parent samples in ﬂow equivalence. Using this bridge function, we predicted the physical properties of the parent samples, including the permeability, bulk porosity, tortuosity, and pore fractal dimension. To validate our model, we performed Lattice Boltzmann (LB) simulations on nine carbonate samples and compared the LB simulation results with our model’s predictions. We also compared our predicted results with available data on various porous materials, such as sandstone, glass beads, and carbonate, in the literature. Our ﬁndings demonstrate that without reconstructing 3D pore networks, our method provides a reliable estimation of sample permeability using 2D cross-sectional images. This approach not only simpliﬁes the determination of sample permeability in heterogeneous porous media but also sheds new light on the inherent correlations between 2D cross-sectional information and 3D pore structures of parent samples. Moreover, the derived model may be conducible to a better understanding of ﬂow in reservoirs during the extraction of unconventional onshore and offshore oil/gas.

Generally, there are three approaches for estimating the permeability of porous material, including experimental measurement, empirical models, and digital rock physics (DRP). Experimental measurement is time-consuming and expensive. Moreover, this approach is not repeatable due to its destructive feature. Instead, empirical models, including the well-known Kozeny-Carman (KC) equation and its variations, have been proposed by various researchers [7,[10][11][12]15,16], of which permeability is related by several static parameters, such as porosity, formation factor, specific surface, and other pore structure parameters. Although these empirical models are efficient and cheap [17], they suffer from general applicability due to the requirement of fitting parameters, which are generally rock-type-dependent [12,18].
To alleviate such limitations, many efforts have been made to directly estimate permeability from available 2D digital images by investigating the similarities between 2D pore structures and corresponding 3D ones [17,25,26,29,30]. Berryman and Blair [31] predicted the permeability of various porous materials by combining the pore information (porosity and two-point spatial correction function) from 2D SEM images with a form of the KC equation and validated their model with experimental data. However, experimental measurements of the formation factor in their model are missing. Yu and Chen [29] derived a fractal model by incorporating pore structure parameters (e.g., fractal dimension and the pore radius) from 2D images. However, direct extraction pore information from 2D images is different from 3D pore information of parent samples, which is supported by Sisavath et al. [32] and Saxena et al. [26]. Thus, transforming 2D pore information into 3D is necessary. Saxena et al. [26] compared rock permeability in 3D and the permeability computed directly from 2D slices and then proposed a bridge function between the 2D results and the 3D permeability. However, there are two fitting parameters (i.e., the normalized magnitude of variation in the pore radius and the "wavelength" of variation in the pore radius) in their model, which are difficult to be obtained. Recently, by combining the KC equation and Young-Laplace equation, Chen et al. [17] derived a bridge function to link the 2D pore structure and 3D pore structure and predicted the capillary pressure curve (CPC) of porous materials. It was found that their approach is useful to predict the CPC of porous media with a 2D pore structure. However, as they suggested, it was of great practical significance to propose other bridge functions to study the corresponding physical properties (e.g., thermal conductivity coefficient, acoustic conductivity) of porous media. To date, how to effectively transform 2D pore information from 2D digital images into 3D pore information without 3D reconstruction and then predict the 3D permeability of porous material still remains an open question and requires further investigation.
In this paper, we developed a novel, robust, and effective approach to predict the 3D permeability of porous materials from 2D images without 3D reconstruction. Previous studies reveal that interspaces in porous media possess fractal characteristics [12,29,[33][34][35][36][37]. Specifically speaking, many scholars [12,29,35,37] suggested that the pore fractal dimension and tortuosity fractal dimension were two key parameters characterizing fluid flow in porous media. Thus, in this paper, these two parameters were applied to characterize the permeability of porous media. As the box-counting algorithm [38] and theoretical models [29,35] can effectively predict the pore fractal dimension [12,[39][40][41], in this paper, the pore fractal dimension was estimated by the box-counting algorithm and theoretical models. In 2015, Wei et al [42]. derived a tortuosity fractal dimension model for fractal porous media that is related to porosity and the pore fractal dimension. As this model is effective at determining the tortuosity fractal dimension [43,44], this model was also used in this paper to model the tortuosity fractal dimension.
Recently, as a promising technique for modeling the fluid flow in porous media at the pore scale, the lattice Boltzmann Method (LBM) has been widely applied by scholars to predict properties (e.g., permeability and relative permeability) of porous media [45][46][47][48][49]. For example, Khodja et al. [46] stated that compared to the traditional methods of computational fluid dynamics, the LBM was ideally suited for massively parallel computation, and it had the advantages of simplicity and flexibility in dealing with complicated geometry. Thus, in this paper, the D3Q19 LBM model was applied to validate our derived model.
In what follows, a new analytical model was derived to transform 2D pore information into a 3D pore structure and give an estimation of the 3D permeability for porous media. Compared with DRP, these models will skip huge workloads of reconstructing 3D pore networks. Subsequently, lattice Boltzmann (LB) simulations were conducted on nine carbonate samples to validate our derived model. In addition, the available experimental data was applied to further validate the developed model and is followed by discussions of this model. Finally, the conclusions are presented.

Methodology
In this section, based on the pore morphology of 2D cross-sectional images, an analytical permeability of 3D global matrix pore structure is presented. It should be noted here that, for simplicity, all the 2D cross-sections are perpendicular to the flow direction; thus, the micro-fractures in actual porous media are ignored. The specific steps are as follows ( Figure 1): Firstly, pore structure of 2D cross-sectional images will be characterized using fractal theory, and the permeability of 2D pore structure K 2D will be derived based on the tube bundle model and fractal geometry [29,39,40]. In the capillary bundle model, the capillary bundles are perpendicular to the images (i.e., capillary bundles run along y direction). Then, by combining the KC equation and fractal geometry, 3D pore structure of porous media in flow equivalence will be estimated, and the permeability of natural samples (in 3D) K 3D will be derived. Finally, a bridge function connecting K 2D and K 3D will be proposed. computational fluid dynamics, the LBM was ideally suited for massively parallel computation, and it had the advantages of simplicity and flexibility in dealing with complicated geometry. Thus, in this paper, the D3Q19 LBM model was applied to validate our derived model.
In what follows, a new analytical model was derived to transform 2D pore information into a 3D pore structure and give an estimation of the 3D permeability for porous media. Compared with DRP, these models will skip huge workloads of reconstructing 3D pore networks. Subsequently, lattice Boltzmann (LB) simulations were conducted on nine carbonate samples to validate our derived model. In addition, the available experimental data was applied to further validate the developed model and is followed by discussions of this model. Finally, the conclusions are presented.

Methodology
In this section, based on the pore morphology of 2D cross-sectional images, an analytical permeability of 3D global matrix pore structure is presented. It should be noted here that, for simplicity, all the 2D cross-sections are perpendicular to the flow direction; thus, the micro-fractures in actual porous media are ignored. The specific steps are as follows ( Figure 1): Firstly, pore structure of 2D cross-sectional images will be characterized using fractal theory, and the permeability of 2D pore structure K2D will be derived based on the tube bundle model and fractal geometry [29,39,40]. In the capillary bundle model, the capillary bundles are perpendicular to the images (i.e., capillary bundles run along y direction). Then, by combining the KC equation and fractal geometry, 3D pore structure of porous media in flow equivalence will be estimated, and the permeability of natural samples (in 3D) K3D will be derived. Finally, a bridge function connecting K2D and K3D will be proposed.

Permeability of 2D cross-section:
In general, based on 2D imaging (thin section, back scatter electron, micro-CT scanning, scanning electron microscope, scanning electron microscopy, and electron probes, etc.) and digital image processing techniques, pore structure parameters (areal porosity, pore size distribution, average pore radius, pore perimeter, specific surface area, pore fractal dimension, coordination number, and shape factor, etc.) of 2D cross-sections can be determined [26,31,50,51]. For instance, pore fractal dimension of 2D cross-section can be obtained by using the box-counting algorithm [38]. In addition, areal porosity of 2D cross-section can be obtained by thin section analysis or digital image processing techniques. Moreover, by using digital image processing techniques, various pore structure parameters (pore size distribution, the maximum pore size, minimum pore size, average pore radius, pore perimeter, and specific surface area) of 2D cross-sectional images can be easily determined [27,36,46,52,53]. For example, for a given 2D digital image, parameters (r max , r min , and ϕ 2D ) can be determined by using the Avizo software, which can effectively determine pore space morphology and extract pore networks [25].
It is well-known that, affected by complex sedimentary and diagenetic processes, pore structure of natural porous media is complicated with obvious heterogeneity, which is difficult to be characterized. Specifically, many scholars suggested that Euclidean mathematics are not suitable to characterize a complex pore structure [29,33,34,54]. With the concept of "Fractal geometry" first proposed by Mandelbrot [33], the problem was well-solved. A large number of studies (experimental and theoretical studies) show that pore structure of most of the sedimentary rocks has obvious fractal phenomena. According to Equation (A4), by combining determined pore structure parameters (e.g., pore radius r, areal porosity ϕ 2D , pore fractal dimension d f ) from 2D imaging and fractal theory, permeability of 2D crosssection (i.e., 2D composite permeated with tortuous capillary tubes that are perpendicular to the 2D cross-section, Figure 1) K 2D is given by [12,29,37,40,55]: where the subscripts max and min denote the maximum and minimum, respectively. More details about Equation (1) can be found in Appendix A. In Equation (1), d t is tortuosity fractal dimension in 2D space, which is [42]: Generally speaking, tortuosity fractal dimension d t ranges from 1 to 2 in 2D space. When the matrix plane is filled with so highly tortuous capillary tubes, d t is equal to 2; however, for straight capillary tubes, d t is equal to 1, and on this condition, Equation (1) can be simplified as: Furthermore, A in Equations (1) and (3) is the cross-sectional area of the representative elementary volume (REV) in 2D space, which can be calculated by Equation (A3) as: Physically, d f in Equation (1) is the pore fractal dimension, which ranges from unity to 2 (1< d f < 2) in 2D pore space. Mathematically, d f can be determined by box-counting algorithm [38] or by the following equation [29,35]: wherein the Euclidean dimension d e is 2 in 2D pore space. Equation (5) reveals that, under a given ratio r min /r max , d f increases with an increasing of areal porosity ϕ 2D . In addition, by fitting the experimental data, Chen et al. [56] also found a positive correlation between d f and ϕ 2D and suggested d f could be estimated by d f = 0.974ϕ 0.173 2D . Equations (1) through (5) form the basis for analysis of K 2D , which were utilized to estimate permeability of 3D parent porous media in flow equivalence.
3D pore structure in flow equivalence: As 2D cross-sections come from 3D parent samples, various studies show that pore structure information (pore morphology) of 2D slice is identical to that of 3D matrix system [17,30,57]. For example, Chen et al. [17] suggested that areal porosity of 2D cross-sections ϕ 2D should be similar to bulk porosity of 3D matrix ϕ 3D . In addition, to estimate permeability in 3D from 2D images, Saxena et al. [26] also assumed that ϕ 2D was similar to ϕ 3D . Figure 2a compares bulk porosity of 3D parent samples versus areal porosity of 2D cross-sections [17,26,31,36,46,53]. Taking the data from Saxena et al. [26], for example, Saxena et al. conducted modeling on 4 samples (i.e., a Fontainebleau sandstone, a Berea sandstone, a Bituminous sand sample, and a Grosmont carbonate). Based on digital rock technology, the bulk porosity ϕ 3D for each sample was determined. Then, multiple 2D cross-sections (lying perpendicular to the flow direction) of a parent cube sample were sampled, and the corresponding areal porosity of each slice was measured using digital image process technology. Similarly, for 3D samples with bulk porosity determined, Wu et al. [36,53] divided the original sample into various layers along vertical direction using low-and high-resolution X-ray CT scanning technologies. Then, areal porosities in each layer of low-and high-resolution digital rock images were measured. Besides digital rock technology, thin section analysis was used by Berryman and Blair [31] and Chen et al. [17] to measure areal porosity of 2D thin. As shown in Figure 2a, porosity of 3D sample ϕ 3D is approximately in the middle between the maximum and minimum values of areal porosities of 2D cross-sections. Specifically, for a 3D porous media with bulk porosity ϕ 3D , most of the areal porosity data ϕ 2D for 2D sectional sections are in the range of 0.9ϕ 3D and 1.1ϕ 3D . Figure 2b compares bulk porosity of 3D samples versus average porosity of 2D cross-sections. Results suggest that porosity of 3D sample is quite consistent with the average areal porosity of 2D slices. Specifically, the slope of the fitting line is 1.009, which is close to unity. And the intercept of the fitting line is close to 3.38 × 10 −4 , which is quite close to 0. As a result, for simplicity of the model, areal porosity of 2D cross-sections ϕ 2D is assumed to be identical to bulk porosity of 3D matrix (parent samples) ϕ 3D , i.e., ϕ 3D = ϕ 2D . In general, for actual pore structure in 3D space, capillary tube is not straight and tortuosity τ (or tortuosity fractal dimension DT) is larger than unity. Based on KC equation, Chen et al. [17] developed a bridge function connecting capillary radius in 2D space to that in 3D space. As stated by Chen et al. [17], for a given capillary with pore radius r In general, for actual pore structure in 3D space, capillary tube is not straight and tortuosity τ (or tortuosity fractal dimension D T ) is larger than unity. Based on KC equation, Chen et al. [17] developed a bridge function connecting capillary radius in 2D space to that in 3D space. As stated by Chen et al. [17], for a given capillary with pore radius r in 2D space, the corresponding pore radius r 3D in 3D space is: As tortuosity τ is larger than unity, Equation (7) reveals that pore radius r 3D in 3D space is less than that in 2D space, and the ratio of r to r 3D is tortuosity to the fourth power, i.e., τ 4 . With the same method as Chen et al. [17], but assuming that the ratio of specific surface in 3D space to that in 2D space is 3/2 and not the tortuosity, Saxena et al. [26] suggested that the ratio of r to r 3D is 9τ 2 /4. To some extent, 9τ 2 /4 can be considered a special case of τ 4 (e.g., when τ is assigned as 3/2, 9τ 2 /4 equals τ 4 ), so, in this paper, to make our model more general, Equation (7) was used to correlate r and r 3D .
According to fractal theory, tortuosity τ of pore radius r 3D is [29]: where V is the bulk volume of the REV. In light of Equation (A10), V can be determined on basis of pore fractal dimension D f in 3D space as: By combining Equations (7) through (9), we have: where the subscripts 3D,max and 3D,min denote the maximum and minimum in 3D space, respectively. In what follows, subscript 3D presents in 3D space. Physically, r 3D,max should be larger than r 3D,min , which means D T in Equation (10) should be less than 1.25. Based on fractal theory, in 3D pore space, D T is [42]: wherein D f is pore fractal dimension in 3D space, which can be also determined by [29,35]: More details about Equation (11) can be found in Equation (A8). By substituting Equation (10) into Equation (12), D f can be rewritten as: By solving Equation (14), parameter D T can be determined. Then, based on Equation (13), D f can be determined. By substituting Equation (10) into Equation (9), we have: By substituting Equation (15) into Equation (10), r 3D,max and r 3D,min can be determined. Then, according to Equation (A11), 3D permeability K 3D of porous media in flow equivalence can be given by: Figure 3 shows the flow chart for the 3D permeability determination. In light of our derived model, the methodology workflow is summarized as follows: Figure 3 shows the flow chart for the 3D permeability determination. In light of our derived model, the methodology workflow is summarized as follows: Step 1: Determine 2D pore structure parameters (e.g., rmax, rmin, φ2D, and df) from 2D digital image. Specifically, for a given 2D digital image, parameters (rmax, rmin, and φ2D) can be determined by using the Avizo software [25]. In addition, df can be determined with box-counting algorithm [38] or with Equation (5). Subsequently, parameter φ3D can be calculated using Equation (6).
Step 2: Determine parameter DT by solving Equation (14). Then, determine Df by using Equation (13). In addition, determine parameter V by using Equation (15). Then, r3D,max and r3D,min can be determined by using Equation (10).
Step 3: Determine K3D and τav of the corresponding parent sample using Equations (16) and (A12). In addition, the bridge function (i.e., the ratio of K2D to K3D) can be determined by combining Equations (1) and (16). Moreover, other 3D pore structure parameters of the parent sample in flow equivalence can be further estimated. Step 1: Determine 2D pore structure parameters (e.g., r max , r min , ϕ 2D , and d f ) from 2D digital image. Specifically, for a given 2D digital image, parameters (r max , r min , and ϕ 2D ) can be determined by using the Avizo software [25]. In addition, d f can be determined with box-counting algorithm [38] or with Equation (5). Subsequently, parameter ϕ 3D can be calculated using Equation (6).
Step 2: Determine parameter D T by solving Equation (14). Then, determine D f by using Equation (13). In addition, determine parameter V by using Equation (15). Then, r 3D,max and r 3D,min can be determined by using Equation (10).
Step 3: Determine K 3D and τ av of the corresponding parent sample using Equations (16) and (A12). In addition, the bridge function (i.e., the ratio of K 2D to K 3D ) can be determined by combining Equations (1) and (16). Moreover, other 3D pore structure parameters of the parent sample in flow equivalence can be further estimated.

Model Validation
Sampling and Experiments: 3 carbonate core plugs (labeled sample 1, sample 2, sample 3) with a diameter of 25 mm and thickness of 6 mm in the Middle East were selected. For these samples, multiple cubical sites were semi-randomly selected and then used in micro-CT X-ray imaging with the Versa XRM-500 X-ray micro-CT scanner and in simulation with the lattice Boltzmann method (D3Q19 model) to determine the permeability in 3D. Taking sample 1, for example, two cubical sites (site A and site B) were semi-randomly selected to avoid vugs. For site A with a bulk porosity ϕ 3D of 23.5% and a total volume of 1.64 3 mm 3 , the volume (voxel) is 648 3 , and the voxel volume is 2.53 3 µm 3 . With the D3Q19 model stated in Appendix B, the determined permeability in 3D is 27.8 × 10 −3 µm 2 . In addition, for site B with a bulk porosity ϕ 3D of 23.4%, the volume (voxel) is 681 3 , the voxel volume is 1.55 3 µm 3 , and the total volume is about 1. To verify the derived model in Section 2, a cross-section of each sample was sliced and used to estimate the permeability of the corresponding parent sample. Specifically speaking, firstly the 2D pore structure parameters (e.g., r max , r min , ϕ 2D , and d f ) from the 2D digital image were determined; then, based on our derived model, the 3D permeability of the carbonates was predicted. In addition, our predicted 3D permeability was compared to that from D3Q19 LBM to validate the feasibility and effectiveness of our derived model. The experimental data and modeling results are summarized in Table 1. As shown in Table 1 and Figure 4a, the areal porosity (2D porosity) of the cross-section consists of the bulk porosity (3D porosity) of the parent sample. The results (Figure 4a) suggest the slope and the intercept of the fitting line are 1.048 (close to unity) and 7.56 × 10 −3 (close to 0), respectively. Similar findings are shown in Figure 2. The results shown in Table 1 suggest that the predicted 3D permeability of each parent sample from the derived model exhibits excellent agreement with the experimental data (the results from the LBM simulations). Taking parent sample 1-A (Figure 4a), for example, in light of the image processing technology, the ϕ 2D of the selected 2D cross-section (Figure 4b) is 31.24%. In addition, with the box-counting algorithm [38], the double logarithmic (Lg) plot of the box size (pixel) versus the total pore number is presented in Figure 4b; then, the pore fractal dimension d f (the negative of the slope) was determined to be 1.686. As shown in Figure 4b, the pore structure of the 2D cross-section of sample 1-A has a typical fractal scale. In our modeling, d f and ϕ 2D were assigned as 1.686 and 31.24%, respectively. r max was assigned as 535.06 µm, and ϕ 3D was assigned as 23.5%. Moreover, r 3D,max and D f were assigned as 41.18 µm and 2.81, respectively. Then, based on the derived model, the predicted permeability K 3D is 27.8 × 10 −3 µm 2 , which is in agreement with that determined by the LBM simulation. Figure 4c presents the pore structure parameters (e.g., τ av , r max /r 3D,max , and D T ) of the samples. As we can see from Figure 4c, the parameter r max is larger than r 3D,max , which is consistent with the results stated by Chen et al. [17]. For these 9 samples, the ratio r max /r 3D,max ranges from 3.18 to 10.58. In addition, D T ranges from 1.10 to 1.15, and τ av ranges from 2.19 to 4.69. The results (Figure 4c) also suggest that even from the same parent sample, the different sites behave strongly in heterogeneity. For instance, sample 1-A and sample 1-B come from the same parent sample 1; however, their physical properties vary greatly. Figure 4d presents the average tortuosity versus the ratio r max /r 3D,max . The results (Figure 4c) show that there exists an approximately linear relationship between the ratio of the maximum pore radius in 2D space to that in 3D space (r max /r 3D,max ) and the average tortuosity of porous media. Physically, when the average tortuosity of porous media is unity, the corresponding ratio r max /r 3D,max approximately equals unity. As we can see from the fitting equation, when τ av is assigned as unity, r max /r 3D,max is determined to be 0.82, which is close to unity, which verifies our model. The reason for the deviation may be caused by the computational error.
In general, carbonates generally contain clay minerals, and the interaction between clay minerals and fluids (e.g., water, gas, oil) will significantly affect the permeability of porous media [37,40]. For example, Lei et al. [37] derived analytical permeability to study the effect of clay swelling on the permeability of clay-rich argillaceous porous media. In addition, as the interaction between solid minerals (e.g., illite, montmorillonite, and kaolinite) and fluids in porous media (e.g., carbonate media) will change pore structures, the effect of the physical or chemical nature of solid minerals on the fluid flow in porous media is of great significance. However, the derived model in this paper does not consider the interaction between solid minerals and fluids. Thus, in our future work, the interaction between solid minerals and fluids will be taken into account to make our model more reasonable.  [26]. For these 12 samples, samples 1 to 3 come from the same parent Berea sandstone and represent Berea sandstones A to C, respectively. In addition, samples 4 to 6 represent Bituminous sands A to C, respectively. Samples 7 to 9 represent Fontainebleau (Font) sandstones A to C, respectively. And samples 10 to 12 represent Grosmont carbonates A to C, respectively. During the numerical experiments, Saxena et al. first calculated the K 3D of the parent samples based on 3D digital rocks using an LBM simulation; then, they sliced 2D thin sections along the flow direction from the parent samples and performed an LBM simulation to estimate the K 2D of the thin section [26]. It should be noted that to calculate the K 2D of the thin section, the pore structure was restructured by permeating the pores in the thin section with straight tubes along the flow direction. That is, for the restructured porous media used to determine K 2D , the tortuosity fractal dimension d t was equal to 1. As a result, during our modeling, d t was assigned as 1, and Equation (3)     Taking sample 1 (sample Berea A) with a volume of 800 × 800 × 800 µm 3 and a pixel size of 0.74 µm, for example, the bulk porosity ϕ 3D is 4%, and the areal porosity ϕ 2D of the thin sections ranges from 2% to 6% with an average of 4%. Based on simulations, the LBM 3D permeability of the parent sample was determined as 3 × 10 −3 µm 2 , and the average LBM 2D permeability of the thin sections is 624 × 10 −3 µm 2 . In our modeling, both ϕ 3D and ϕ 2D were assigned as 4%. Based on Equations (3) and (13), K 2D and K 3D can be calculated. The results (Figure 5a-d) suggest that the predictions (K 2D and K 3D ) provide a good match over the test data. Besides K 2D and K 3D , other parameters (e.g., d f , τ av , r max /r 3D,max , D T , and D f ) of each sample are presented in Figure 5e. As we can see from Figure 5e for these 12 samples, τ av ranges from 1.10 to 1.50, which is smaller than that of the carbonates in Figure 4. This reveals that the carbonates in Figure 4 have stronger heterogeneity than the rocks in Figure 5. The reason may be that the rocks in Figure 4 are carbonate core plugs; however, most of the rocks in Figure 5 are sandstones, which are more homogeneous. In addition, Figure 5e shows that the ratio r max /r 3D,max ranges from 1.15 to 1.74, D T ranges from 1.01 to 1.03, d f ranges from 1.64 to 1.82, and D f ranges from 2.69 to 2.83. The results (Figure 5f) also suggest that τ av and the ratio r max /r 3D,max are linearly dependent. When τ av is assigned as unity, the ratio r max /r 3D,max is estimated to be 0.93, which is approximately equal to unity. Similar findings are shown in Figure 4e.
Dataset of Berryman and Blair [31]: Based on digital image analysis, Berryman and Blair first studied the two-point correlation functions, porosity, and specific surface area of glass beads (55 µm), Ironton Galesville (IG) sandstones, and Berea sandstones [31]. Then, they predicted the permeability of these samples with the following equation: where F is the formation factor, and b is a constant that depends on the cross-section of the tubes. It is recommended that b is assigned as 2 and 3 for circular tubes and flat cracks, respectively. In order to further verify our derived model, the predictions from our derived model were compared against the available data of Berryman and Blair [31]. Taking sample 3 (sample Berea 1) with an image specific surface of 0.0281 µm −1 , an areal porosity ϕ 2D of 17%, and an image permeability K 2D of 0.312 µm 2 . For example, Berryman and Blair found the bulk porosity (laboratory porosity) ϕ 3D ranged from 15% to 18% with an average value of 16.5%. In addition, based on Equation (14), the laboratory permeability of K 3D was determined to be 0.023 µm 2 . In our modeling, ϕ 3D and ϕ 2D were assigned as 16.5% and 17%, respectively. The results (Figure 6a) suggest that the predicted results (K 2D and K 3D ) of these samples (i.e., two glass beads, four Berea samples, IG-775 and IG-785) provide a good match over the test data. In addition, the pore structure parameters (e.g., d f , τ av , r max /r 3D,max , D T , and D f ) of each sample are presented in Figure 6b. As pore structures (d f , τ av , r max /r 3D,max , D T , and D f ) of these samples are different, the curves show different behaviors. In addition, the results (Figure 6c) also suggest that τ av is remarkably correlated linearly with the ratio r max /r 3D,max .
For a given porous medium (bulk porosity ϕ 3D ) composed of identical spherical grains, the average particle radius R ap of the porous medium could be written as [26,31,58]: where S 3D is the specific surface area. Based on the derivation in Appendix C, R ap can be also expressed as: where λ is the weight coefficient of the average spherical particle radius, R avs is the average spherical particle radius, and R avc is the average circular particle radius. More details about R avs and R avc can be found in Equations (A22) and (A24). By combining Equations (15) and (16), the weight coefficient of the average spherical particle radius can be determined. Figure 6d compares the calculated average particle radius for spherical particles and the calculated average particle radius for circular particles with that from the former model. The results from Figure 6d demonstrate that the calculated average particle radius from Equation (15) is approximately in the middle between the calculated values from Equations (A22) and (A24). This indicates that the model gives predicted values that are quite consistent with the results from the former model. In addition, by combining the results from Equations (15) and (16), the weight coefficients for different samples have been determined in Figure 6d. Taking sample 3 for example, when the weight coefficient is unity, the average particle radius can be effectively determined by Equation (A22).    For a given porous medium (bulk porosity φ3D) composed of identical spherical grains, the average particle radius Rap of the porous medium could be written as [26,31,58]: (18) where S3D is the specific surface area. Based on the derivation in Appendix C, Rap can be also expressed as: where λ is the weight coefficient of the average spherical particle radius, Ravs is the average spherical particle radius, and Ravc is the average circular particle radius. More details about Ravs and Ravc can be found in Equations (A22) and (A24). By combining Equations (15)   (c) τ av versus the ratio r max /r 3D,max for various samples; (d) the predicted average particle radius for spherical particles and circular particles versus that from Equation (15).

Results and Discussions
After validation with exhaustive experimental data, this derived model was utilized for sensitivity analysis of different parameters (D f and D T ) on the 3D permeability of porous media. Figure 7 shows the influence of pore fractal dimension D f on the average tortuosity τ av and K 3D of porous media in 3D. In this case, the maximum and minimum pore radii in 2D were 3 µm and 0.13 µm, respectively. The parameters ϕ 2D and ϕ 3D were both 0.15, and parameter D T was 1.1. During the calculation, the parameter D f ranged from 2.2 to 2.8. Based on our derived model, the maximum pore radius in 3D space r 3D,max was determined, which ranged from 1.4 µm to 0.82 µm. In addition, the determined parameter r 3D,min ranged from 0.0077 µm to 0.0045 µm. Specifically, when parameter D f was assigned as 2.2, the corresponding values of r 3D,max and r 3D,min were 1.4 µm and 0.0077 µm, respectively. As one can see from Figure 7a, there exists a positive relationship between D f and the average tortuosity τ av . The main reason is that a larger value of D f means a more complex pore structure of porous materials, leading to a larger value of τ av . Furthermore, the permeability of porous media in 3D decreases with an increase in D f , which is anticipated. Figure 7b demonstrates that the linear correlativity between τ av and the ratio r max /r 3D,max is very prominent. Moreover, when τ av is assigned as unity, r max /r 3D,max is determined to be 0.9951, which is close to unity. determined, which ranged from 1.4 µ m to 0.82 µ m. In addition, the determined parameter r3D,min ranged from 0.0077 µ m to 0.0045 µ m. Specifically, when parameter Df was assigned as 2.2, the corresponding values of r3D,max and r3D,min were 1.4 µ m and 0.0077 µ m, respectively. As one can see from Figure 7a, there exists a positive relationship between Df and the average tortuosity τav. The main reason is that a larger value of Df means a more complex pore structure of porous materials, leading to a larger value of τav. Furthermore, the permeability of porous media in 3D decreases with an increase in Df, which is anticipated. Figure 7b demonstrates that the linear correlativity between τav and the ratio rmax/r3D,max is very prominent. Moreover, when τav is assigned as unity, rmax/r3D,max is determined to be 0.9951, which is close to unity. The influence of parameter DT on the τav and K3D are shown in Figure 8. In this case, the basic parameters applied in the model are summarized in the corresponding figures. During the calculation, the parameter Df was assigned as 2.2, and the parameter DT ranged from 1.05 to 1.2. Based on our derived model, the parameter r3D,max ranges from 2.2 µ m to 0.07 µ m. The results (Figure 8a) suggest that a larger τav corresponds to a larger DT. However, K3D decreases as DT increases. Specifically, for this case, when DT increases up to a certain value (e.g., DT ≥ 1.2), the value of K3D is extremely small. The main reason is that a  The influence of parameter D T on the τ av and K 3D are shown in Figure 8. In this case, the basic parameters applied in the model are summarized in the corresponding figures. During the calculation, the parameter D f was assigned as 2.2, and the parameter D T ranged from 1.05 to 1.2. Based on our derived model, the parameter r 3D,max ranges from 2.2 µm to 0.07 µm. The results (Figure 8a) suggest that a larger τ av corresponds to a larger D T . However, K 3D decreases as D T increases. Specifically, for this case, when D T increases up to a certain value (e.g., D T ≥ 1.2), the value of K 3D is extremely small. The main reason is that a larger value of D T means larger seepage resistance, resulting in a small value of K 3D . Figure 8b also reveals that τ av has distinct linear correlations with the ratio r max /r 3D,max . Similar findings have been also demonstrated in Figure 6.  Figure  8b also reveals that τav has distinct linear correlations with the ratio rmax/r3D,max. Similar findings have been also demonstrated in Figure 6. Based on the definition of the bridge function fbridge connecting K2D and K3D (i.e., fbridge = K2D/K3D), we studied the influences of Df and DT on this bridge function fbridge. Figure 9 presents the influences of Df and DT on fbridge. For the calculation, the parameters applied in the model are summarized in the corresponding figures, which are identical to those for Figures 7 and 8, respectively. As one can see from Figure 9, the bridge function increases as Df (or DT) increases. The main reason is that for a given 2D pore structure, a larger value of Df (or DT) corresponds to a smaller value of K3D. Similar findings can be found in Figures 7a and 8a. Moreover, Figure 9 reveals that with an increase of Df (or DT), the increase rate of fbridge increases. The main reason is that with an increase of Df (or DT), the increase rate of the seepage resistance in the porous media increases. Based on the definition of the bridge function f bridge connecting K 2D and K 3D (i.e., f bridge = K 2D /K 3D ), we studied the influences of D f and D T on this bridge function f bridge . Figure 9 presents the influences of D f and D T on f bridge . For the calculation, the parameters applied in the model are summarized in the corresponding figures, which are identical to those for Figures 7 and 8, respectively. As one can see from Figure 9, the bridge function increases as D f (or D T ) increases. The main reason is that for a given 2D pore structure, a larger value of D f (or D T ) corresponds to a smaller value of K 3D . Similar findings can be found in Figures 7a and 8a. Moreover, Figure 9 reveals that with an increase of D f (or D T ), the increase rate of f bridge increases. The main reason is that with an increase of D f (or D T ), the increase rate of the seepage resistance in the porous media increases.
Based on the definition of the bridge function fbridge connecting K2D and K3D (i.e., fbridge = K2D/K3D), we studied the influences of Df and DT on this bridge function fbridge. Figure 9 presents the influences of Df and DT on fbridge. For the calculation, the parameters applied in the model are summarized in the corresponding figures, which are identical to those for Figures 7 and 8, respectively. As one can see from Figure 9, the bridge function increases as Df (or DT) increases. The main reason is that for a given 2D pore structure, a larger value of Df (or DT) corresponds to a smaller value of K3D. Similar findings can be found in Figures 7a and 8a. Moreover, Figure 9 reveals that with an increase of Df (or DT), the increase rate of fbridge increases. The main reason is that with an increase of Df (or DT), the increase rate of the seepage resistance in the porous media increases.

Model advantages and limitations:
The developed model provides a theoretical basis for predicting the 3D permeability of porous media from 2D digital image analysis without the need for 3D reconstruction, achieving high accuracy even with only 2D information. Compared to DRP, our model not only reduces the computational cost of high-resolution scanning, 3D pore network reconstruction, and numerical simulation but also captures the effect of realistic pore shapes on permeability. Furthermore, our model can be used for inverse modeling to estimate relevant parameters, such as tortuosity and the pore fractal dimension in 3D, making it highly practical for estimating permeability in heterogeneous porous media. Overall, our derived model is of great practical significance, as 2D images are easier and cheaper to access than 3D images and reconstructions of 3D pore networks.
However, it is worth noting that our derived model shares the problems of uncertainty and challenge in terms of predicting 3D pore structures and the permeability of porous media. Physically speaking, for a given 2D image, the possible 3D samples are infinite. That is, for a given 2D image, there are infinitely many permeabilities to this problem. As our model is derived based on fractal theory, it is limited to fractal porous media; however, it may not be suitable for some porous media. In addition, in our modeling, the areal porosity of 2D cross-sections ϕ 2D is assumed to be identical to the bulk porosity of the 3D matrix. However, ϕ 2D varies with the slices, and sometimes the difference between the maximum and the minimum 2D porosity is relatively large. Although ϕ 3D is approximately in the middle between the maximum and minimum values of ϕ 2D , 2D porosity is not enough to represent ϕ 3D . Moreover, our developed model is limited to intact porous media and ignores the effect of micro-fractures on the 3D permeability of porous media. Thus, further research is required to reduce the uncertainty in estimating the 3D permeability of porous media from 2D images without reconstruction. Furthermore, in general, the pore surface of porous materials is rough, and the surface fractal dimension is crucial to characterize the fluid flow in porous media [41,59,60]. For example, Lei et al. [59] and Xiao et al. [41] derived theoretical models to study the fluid flow in porous media, and they concluded that the effects of surface fractal dimension on permeability and relative permeability in porous media were significant. Thus, to make our model more reasonable, in our future work, a rough pore surface and surface fractal dimensions will be taken into account. Moreover, as mentioned above, the interaction between solid minerals and fluids will significantly affect the permeability of porous media. To improve the applicability of our derived model, the interaction between solid minerals and fluids will be taken into account. What is more, as the information on fracture systems may not be obtained from 2D cross-sections, in this paper, our derived model focuses on predicting the 3D permeability of intact porous media, and the fractures are ignored. In our future work, we will try to extend our model to study the 3D permeability of fractured porous media from 2D cross-sections of parent samples without 3D reconstruction.

Conclusions
In this paper, a novel analytical model was derived to estimate the 3D permeability of porous media from 2D cross-sectional images without reconstruction. Based on fractal theory and the Kozeny-Carman equation, a bridge function was developed to correlate the 2D pore information from the 2D images and the 3D pore structure of the parent samples. The derived model was validated with the results from lattice Boltzmann method (LBM) simulations and various experimental data and is shown to perform well.
The derived model was also used to conduct sensitivity analysis of different parameters (e.g., tortuosity fractal dimension in 3D space D T , pore fractal dimension in 3D space D f ) on the 3D permeability K 3D of porous media. The results indicate that the average tortuosity τ av decreases with an increase of D f (or D T ). In addition, the average tortuosity τ av is remarkably correlated linearly with the ratio of r max (the maximum pore radius of 2D images) to r 3D,max (the maximum pore radius of 3D porous media). Moreover, the permeability K 3D decreases as D f (or D T ) increases. With an increase of D f (or D T ), the decrease rate of the permeability K 3D increases.
The proposed model not only reveals the Intrinsic link of 2D pore information and 3D pore structure in unconventional oil/gas reservoirs but also reduces the computational cost of high-resolution scanning and flow simulation to predict the 3D permeability of unconventional reservoirs. Our proposed model can be applied in many fields, such as CO 2 geology storage, unconventional onshore and offshore oil/gas development, and groundwater seepage.
Although the model focuses on predicting single-phase flow in dry porous media, when the irreducible wetting phase saturation is taken into account, our model can be extended to study the 3D effective permeability of porous media. In addition, it is also available to extend this model to study multi-phase flow in porous media and obtain relative permeability in porous media.
However, it should be noted that our proposed model never considers a rough pore surface (e.g., surface fractal dimension) or interactions in solid minerals-fluid systems. Thus, in our future work, more mechanisms will be taken into account to make our model more reasonable.

Conflicts of Interest:
We declare that we have no conflict of interest.

Latin symbols A
The cross-sectional area of the representative element (µm 2 ) A av The average pore area (µm 2 ) b A constant that depends on the cross-section of the tubes (dimensionless) c Lattice sound speed, which is determined by ∆x/ ∆t (m/s) d e The Euclidean dimension is 2 in 2D space (dimensionless) d f The pore fractal dimension in 2D space (dimensionless) d t The tortuosity fractal dimension in 2D space (dimensionless) D f The pore fractal dimension in 3D space (dimensionless) D T The tortuosity fractal dimension in 3D space (dimensionless) e α Lattice velocities (m/s) e y The unit vector along the y-axis (m/s) F The formation factor (dimensionless) f The pore size distribution in 2D space (dimensionless) f 3D The pore size distribution in 3D space (dimensionless) f α (y, t) The evolution of the density distribution function (kg/m 3 ) K 2D Permeability of the 2D cross-section (10 −3 µm 2 ) K 3D 3D permeability of porous media in flow equivalence (10 −3 µm 2 ) K ij Permeability in the i direction when the flow is driven in the j direction (10 −3 µm 2 ) L The actual streamlined length of a tortuous capillary (µm) L j The size of the computational domain in the j direction (µm) L u The edge length of the cubical unit cell (µm)

N
The total number of pores in the representative element of a 2D cross-section (dimensionless) N 3D The total number of pores in a representative elementary volume (dimensionless) N total The total lattice number (dimensionless) p The transient pressure (Pa) p in The inlet transient pressure (Pa) p out The outlet transient pressure (Pa) r The pore radius in 2D space (µm) r av The average pore radius in 2D space (µm) r 3D The pore radius in 3D space (µm) r 3D,av The average pore radius in 3D space (µm) r max The maximum pore radius in 2D space (µm) r min The minimum pore radius in 2D space (µm) r 3D,max The maximum pore radius in 3D space (µm) r 3D,min The minimum pore radius in 3D space (µm) R The particle radius of the porous media (µm) R ap The average particle radius of the porous media (µm) R avc The average circular particle radius (µm) R avs The average spherical particle radius (µm) R max The maximum particle radius (µm) S 2D The specific surface area of the 2D cross-section (µm −1 ) S 3D The specific surface area of the 3D cross-section (µm −1 ) u j The velocity at the void point j (m/s) u The flow velocity (m/s) Then, by the definition of specific surface area S 3D in 3D space, we have: where u i is the component in the i direction of the volumetric average velocity and computed by a summation of u i over the void/fluid grids, N total is the total number of both void and solid grids, K ij is the component in the i direction when the flow is driven in the j direction, and L j is the size of the computational domain in the j direction. The permeability K ij as a tensor has nine components in total by driving the flow in three different directions using three independent simulations. However, the current simulations in this paper always drive the flow by a pressure difference along the j = y direction and focus on the component in the i = y direction, and K yy is denoted by K 3D as the 3D permeability of the porous media. It is worth noting that attention has been paid to alleviating the velocity and permeability dependences on the relaxation time due to the nonzero Knudsen effect in the LBM simulations [46,68] where the Knudsen (Kn) number is √ π/6(τ 0 − 0.5)/N pore [69], and N pore is the grid number used to discretize the dominant pore size. However, we could not set (τ 0 − 0.5) → 0 for Kn → 0 ; that would result in zero kinematic viscosity and an infinite Reynold (Re) number, leading to the dependence of the permeability on Re. In practical simulations of digital rocks [46,68], using (τ 0 − 0.5) around 0.2 and a small relative pressure drop of only 1% (i.e., p in − p out = 0.01p in ) can make the Kn and Re effects negligible, respectively. On the other hand, Li et al. (2018) found the multi-relaxation time (MRT) model might provide a parameter range that was wider than that of the current BGK model, but the basic mechanism remained the same; i.e., in simulations of intrinsic permeability, one should choose the appropriate relaxation parameter to ensure that Kn is small (but also nonzero as in the current BGK model).

Appendix C. Determination of the Average Particle Size
In this section, the average rock particle radius in porous media with bulk porosity ϕ 3D was estimated with the idealized geometrical model shown in Figure A1, assuming that the unit cell was composed of spherical particles of the same size. Based on Figure A1, the volume of the cubical unit cell V u and its edge length L u are: permeability, one should choose the appropriate relaxation parameter to ensure that Kn is small (but also nonzero as in the current BGK model).

Appendix C. Determination of the Average Particle Size
In this section, the average rock particle radius in porous media with bulk porosity φ3D was estimated with the idealized geometrical model shown in Figure A1, assuming that the unit cell was composed of spherical particles of the same size. Based on Figure  A1, the volume of the cubical unit cell Vu and its edge length Lu are: Based on Equation (A18), the average pore volume Vav is: Mathematically, the average pore volume Vav can be also estimated with: Based on Equation (A18), the average pore volume V av is: Mathematically, the average pore volume V av can be also estimated with: wherein r 3D,av is the average pore radius in 3D space, which is: By combining Equations (A19) and (A20), the average spherical particle radius R avs is By simplifying spherical particles in Figure A1 as circular particles, Xu and Yu derived an analytical model for the maximum particle radius R max , which is [12]: In light of Equation (A23), the average circular particle radius R avc can be also determined as: Then, by using the average particle radius statistical average method, the average particle radius R ap of the porous media can be expressed as: wherein λ is the weight coefficient of the average spherical particle radius.