Reconstruction of 3D multi-mineral shale digital rock from a 2D image based on multi-point statistics

Introduction: Shale oil and gas reservoirs contain a variety of inorganic and organic pores that differ significantly from conventional reservoirs, making traditional experiments ineffective. Instead, the pore-scale imaging and modeling method, regarded as a novel and practical approach, is proposed to characterize shale microstructure and petrophysical properties. Therefore, it is of great significance to accurately reconstruct the three-dimensional (3D) microstructure of the porous medium, that is, the digital rock. However, microstructural images of shale at high-resolution, obtained through scanning electron microscopy (SEM) are constrained in the two-dimensional (2D) scale. Method: In this work, a novel iterative algorithm to reconstruct 3D multi-phase shale digital rock from a 2D image using multi-point statistics has been proposed. A multi-grid data template was used to capture the conditional probabilities and data events. The novelty of this work stems from an accurate representation of different types of pores and the mineral characteristics of shale rock from 2D images. Result: A series of simulations were conducted to reconstruct 2D shale digital rock from a 2D segmented training image, 3D shale digital rock from a 2D segmented training image, a 2D gray training image to reconstruct 2D shale digital rock, and a 2D gray training image to reconstruct 3D shale digital rock. Discussion: To corroborate the accuracy of the reconstructed digital rock and evaluate the reliability of the proposed algorithm, we compared the construction image with the training image with the two-point correlation function, geometry, morphological topology structure, and flow characteristics. The reconstruction accuracy indicates that the proposed algorithm can replicate the higher-order statistical information of the training image.


Introduction
Shale reservoirs are distributed worldwide and are gaining more attention (Ma et al., 2018;Feng et al., 2020;Kang et al., 2020;Solarin et al., 2020;Katz et al., 2021;Jiang et al., 2022). The microstructure of shale is complex and includes different types of inorganic pores, porous and non-porous organic matter, and various mineral phases. Furthermore, different phases and their shapes are distributed randomly and disorderly. Moreover, the pore structures of shale range from nanometers to micrometers and exhibit multi-scale features, which make the petrophysical properties and flow mechanism of shale thoroughly distinct compared with conventional reservoirs (Josh et al., 2012;Yang et al., 2013;Sakhaee-Pour and Bryant, 2015;Deng et al., 2016;Ilgen et al., 2017;Huang et al., 2018). Therefore, traditional, widely used experimental methods are difficult to apply directly to shale research.
In the last decade, with the increase in computing power and the development of advanced imaging equipment, pore-scale imaging and modeling combined with numerical algorithms have become a new promising and outstanding petrophysical technology for imaging the microstructure of porous media, which is utilized to calculate the rock properties such as porosity, permeability, and relative permeability (Andrä et al., 2013a;Andrä et al., 2013b). Because of their advantages in terms of speed, accuracy, low cost, and reuse, pore-scale imaging and modeling have been widely used in petroleum engineering (Blunt et al., 2013). This process involves obtaining the shale rock microstructure and modeling transport in porous media (Bultreys et al., 2015). It is very important, while building 3D digital rocks, to get an accurate 3D microstructure of porous media. This is because the reconstructed porous structure is closely related to petrophysical properties and flow characteristics. At present, two main methods can be applied for digital rock reconstruction, namely physical experimental and numerical reconstruction approaches (Zhu et al., 2019). Advanced optical instruments were used in the physical experiments to obtain the microstructure of the rock samples. These include confocal laser scanning (Hackley and Kus, 2015;Hackley et al., 2020), X-ray micro-and nano-computed tomography (XCT) scanning (Gu et al., 2015;Tang et al., 2016;Zhou et al., 2016), and focused ion beam scanning electron microscopy (FIB-SEM) imaging (Cnudde and Boone, 2013;Ma et al., 2017;Zhang et al., 2020). Although physical experiments can accurately obtain the actual microstructure of a rock sample, the method is timeconsuming and expensive and cannot balance image resolution with rock sample size. Due to these disadvantages, it is difficult to establish a 3D model for heterogeneous shale rock with high resolution using physical experiments. SEM is easily and widely used to reveal the microstructure of shale in a two-dimensional (2D) format with high resolution (Klaver et al., 2016;Saif et al., 2017). As a consequence, it is of great significance to develop accurate, stable, and effective 3D digital rock reconstruction algorithms based on 2D slices.
Numerical reconstruction employs a small amount of information to reconstruct 3D digital rocks using mathematical methods. In general, they can be divided into three types. The first is referred to as a "process-based reconstruction method," which simulates the real settlement, compaction, and diagenesis processes of rocks (Tian et al., 2019). It has the advantage of better connectivity in reconstructed 3D digital rocks, and the pore structure is ideal. This is not suitable for shale, which has a complex diagenetic environment. The second method includes the Gaussian field simulation method (Hyman and Winter 2014), the simulated annealing method (SA) (Gerke et al., 2012;Shojaeefard et al., 2016), and Markov Chain Monte Carlo (MCMC) (Wu et al., 2006), which reconstructs 3D digital rocks with statistical information such as the connectivity function obtained from 2D images. The 3D digital rocks constructed according to the Gaussian field simulation approach exhibit poor connectivity. The SA method is suitable for reconstructing 3D digital rocks with large porosity, high permeability, and a simple microstructure. Both computational efficiency and accuracy decrease with an increase in microstructure complexity. The MCMC approach is feasible for building reliable 3D digital rocks. However, it failed to reproduce the microstructure of the anisotropic porous medium. These two methods use low-order statistical information from training images. Thus, it is difficult to reconstruct the global connectivity of the porous medium. The last method for establishing 3D digital rocks is to use high-order characteristic information obtained from training images as well as multi-point statistics (MPS) with spatial uncertainty interpolation (Wu et al., 2018). Compared to the first two methods, the MPS method uses high-order statistical information that reflects the microstructure of the training images. Thus, the adoption of the MPS method to construct the 3D digital rock can be viewed as a feasible approach to achieving the best long-term connectivity. Additionally, newer machine learning algorithms (Zha et al., 2020;Wang et al., 2021;Shan et al., 2022) for reconstruction are also widely used, and convolutional deep belief networks (CDBN) are used to construct stochastic models of porous media (Cang et al., 2017). To reconstruct relatively homogeneous sandstone and heterogeneous carbonate samples, generative adversarial networks and related variant-based methods are used (Mosser et al., 2017;Feng et al., 2019;Shams et al., 2020). Algorithms require a large amount of data for training, consume a lot of memory, limit the size of the reconstructed volume, must be retrained for different types of rocks, and cannot be migrated.
The MPS method for reconstructing 3D digital rocks includes three steps. First, we need to define training images and data templates. Second, it is to sample with the defined training images and data templates. The last step is to recover the microstructure based on Monte Carlo simulations and interpolation. One can divide sampling and recovery methods into two groups: those that use pixels and those that use patterns. Okabe and Blunt first proposed the MPS approach for establishing 3D digital rocks from 2D training images that are isotropic in each direction (Okabe and Blunt, 2005). Because the MPS method is restricted by its poor reconstruction efficiency, many scholars have transplanted MPS reconstruction from CPU devices to GPU devices or developed other algorithms to speed it up (Huang et al., 2013;Zhang et al., 2015). An excellent MPS approach on the basis of the cross-correlation function geostatistical method (CCSIM) to reconstruct 3D digital rock was put forward by Tahmasebi and Sahimi, (2012). The isomapbased MPS that adopts the isomap approach for reducing dimensions and redundancy was developed by Zhang et al, (2016). Also, they came up with a pattern-based method called FILTERSIM. This method uses filters to reduce the number of dimensions in patterns and put them into finite classes (Zhang et al., 2006). In practical applications, it is often difficult to reconstruct 3D porous media using only 2D training images. One solution is to build 3D porous media by rotating 2D training images and establishing probabilities in different directions to achieve probability fusion, but this assumes that all directions are isotropic and cannot solve the problem of anisotropy (Comunian et al., 2012). An alternate solution is to adopt the MPS algorithm to construct systematic images based on 2D training images to form a 3D porous medium by stacking, which leads to the final porous medium reproducing noise and discontinuous structures . To integrate the advantages of both methods, a hybrid approach generates multiple 2D slices as conditional data and 3D training images (Tahmasebi et al., 2015), and the rest of the region is simulated by the MPS algorithm. However, it is challenging to select an applicable switching point between these two approaches, and there are still drawbacks. In this study, an iterative algorithm has been proposed that can use a single training image to rebuild 3D digital rocks.
Even though the MPS method is the most common way to reconstruct 3D two-phase (pore and grain phase) digital rock for conventional reservoirs like sandstone, it still needs to be improved for unconventional reservoirs like shale reservoirs, especially for reconstructing 3D multi-phase digital rocks (Yang et al., 2015;Liu et al., 2018;Yan et al., 2018). Nie et al. came up with a multi-phase reconstruction method with the MCMC reconstruction algorithm to obtain a large-scale multi-phase 3D digital rock (Nie et al., 2016). They also adopted the reconstructed 3D digital rock to perform electrical simulations. Zhu and Yu proposed an approach to incorporate cemented phases into digital sandstone rocks that have already been reconstructed (Zhu et al., 2012). Based on digital rock physics, Saad et al. could reconstruct a multi-phase porous rock (Saad et al., 2018). Most methods proposed to generate multi-phase digital rocks are reconstructed  and then nested directly, without correlating the geology and minerals (for instance, organics and pyrite are frequently found inn correlation with one another). In conclusion, it is a great challenge to reconstruct multi-phase 3D digital rocks, and more research needs to be conducted in the future. Given the various types of pores and mineral characteristics found in shale, developing multi-phase digital rocks from shale is critical and must be pursued further.
In this study, we adopted multi-point statistics to generate a 3D multi-phase shale digital rock considering minerals, different types of pores, and organic components from 2D images. The remaining parts of this paper are structured as follows: In Section 2, the MPS method and workflow for reconstructing the 3D multi-mineral shale digital rocks are explained. Section 3 introduces the evaluation methods for the reconstruction of shale digital rocks. Section 4 depicts various reconstruction cases, such as a 2D segmented training image to reconstructed 2D shale digital rocks, a 2D segmented training image to reconstructed 3D shale digital rocks, and a 2D grayscale training image to reconstructed 3D shale digital rocks. The study is summarized in Section 5.

Methodology of multi-point statistics
In this section, the basic theory related to MPS is introduced. The basic principle for reconstructing digital rocks based on MPS is briefly described to explain why the simulation spaces can Illustration of (A) an image used for training purposes, split into a grain phase (black) and a pore phase (white), (B) a 9×9 2D data template, and (C) a data event acquired after a 9×9 2D data template scanned the training image.

Frontiers in Earth Science
frontiersin.org reproduce the microstructure of the training images. A reconstruction workflow is also presented.

Basic concepts of multi-point statistics
MPS were proposed using traditional two-point statistics based on variograms. The key components include training images, data templates, and data events (Caers and Zhang, 2004;Pyrcz and Deutsch, 2014). The training images contain various feature patterns to be reconstructed in the simulations, as shown in Figure 1A). The black region indicates the grain phase, while the white one depicts the pore phase. It is a collection of conceptual feature patterns that require stability instead of high precision or specific distribution. The training image directly determines the reconstruction quality.
A data template consists of n vectors, which is set to τ n ={h α ; α=1, 2, . . . , n}. The point u in the template is the center, while the additional positions are denoted by the notation u α =u+ h α , where α can represent any value from one to n. Figure 1B shows a 2D data template consisting of 9×9 nodes, where u α is decided by the center point u and 80 vectors h α .
For a data template with attribute S representing the pore, grain, mineral or organic phases, it can take m state values {s k ; k = 1, 2, . . . , m}. The "data event" d(u) consisting of n state values of the n vector u α positions in the data template, is expressed as: where S (u α ) represents the state value at the u α position, d(u) represents the S (u 1 ), . . ., S (u n ) of the n vectors, and is the state value s k1 , . . ., s kn . Figure 1C) shows one of data events captured by the two dimensional data template of Figure 1B) to scan the training image of Figure 1A).

Basic principles for reconstructing digital rock based on multi-point statistics
Data templates are used to scan training images to capture the data events, which can be considered as microstructure features of the training data. It is possible to obtain all of the data events contained in the training image if the training images are scanned while simultaneously moving data templates from the top left to the bottom right node, one node at a time. During the course of scanning the training image with any data template, a data event (composed of the simulated point data u and the condition data in its data template) is captured by the data template. As to the training image, if a data event is repeated N times, the number of repetitions is N. For any point u to be simulated, it is necessary to decide on a function referred to as the conditional probability distribution function (CPDF) of attribute S(u) which can be expressed by the Bayesian conditional probability formula, for any K state value at the condition of a given n data value S (u α ) as: where the denominator is the pattern occurrence probability, and the numerator is the probability that a pattern and a state value to be simulated occur simultaneously. Therefore, the above expression can be rewritten as: Here, c (d n ) stands for the total number of times a particular pattern is repeated, and c k (d n ) is the number of points S(u) to be simulated in the repeated pattern, equal to s k . The "search tree" consists of a set of nodes, where each node in the tree data structure is a data event. The memory space and the scanned training images can be reduced when the CPDF are stored using the data structure of the "search tree." Using numerous data templates is a viable option for extracting detailed information about the pore structure at multiple scales from the training image. This method not only improves sampling efficiency but also allows for the collection of large-scale porous geometry from training images. The MPS method has mainly been applied in the field of macroscopic geological modeling since it was proposed. For digital rock reconstruction, the MPS method is applied to homogeneous and two-phase (pore and grain) sandstone rocks, based mainly on the pixel method. This approach is ineffective for reconstructing the shale digital rock, particularly multi-scale and multi-phase shale digital rock. The microstructure of shale is heterogeneous, and the long-distance connectivity functions of the different phases are inconsistent. In the process of reconstructing a digital rock using the MPS technique, it may be necessary for recovering the pore structure of the training images via spatial interpolation in the area to be simulated. A pixel-based method is used to recover the pore structure of the training images. Based on the microstructure information of the training images, the areas were simulated using the Monte Carlo method, restricted by the CPDF and data events. Once the area to be simulated was traversed, the spatial structure of the entire area to be simulated was reconstructed. Because the probability estimation method was adopted, the simulation results were random. In addition, the entire process was based on the information from the training images, so it could reveal the various possible distributions of the state values in space very well.
2.3 Procedures for reconstructing shale 3D multi-phase digital rock using the multi-point statistics method Step 1: Input a representative 2D shale image, which is a grayscale image or segmented image with different phases. The Frontiers in Earth Science frontiersin.org 04 training image can be obtained using different scanning equipment such as XCT, SEM, etc.
Step 2: Set the representative 2D shale image as the training image. Multiple data templates were used to scan and store the obtained data events and CPDF in the search tree.
Step 3: Specify a random path to access the path of the simulated area. A selected phase (such as the pore, grain, or another phase) is randomly placed in the path, which is treated as the starting point of the simulation. The CPDF is updated and the next simulation point based on a random path using the Monte Carlo method is generated.
Step 4: Determine the dimensionality of the space to be simulated, as shown in the workflow in Figure 2. If the simulation space is 2D, the 2D digital rock will be generated directly based on steps 1-3.
Step 5: If the simulation is 3D, based on Step 4, a 2D selection will be generated and regarded as a new training image, as shown in the workflow in Figure 2. Repeat steps one to four to generate another new cross-section. The 3D digital rock was reconstructed from the series of 2D cross-sections one by one on the basis of the interactive approach.

Evaluation method
After reconstructing the multi-mineral shale digital rock, it is necessary to analyze its quality by evaluating

FIGURE 2
Flow chart of multi-mineral shale digital rock reconstruction using MPS.

Frontiers in Earth Science
frontiersin.org the divergence between the reconstructed results and the training image of the shale. The feasibility of this reconstruction method should also be evaluated. Then, the two-point correlation function, the petrophysical properties, and the pore structure would be compared. The MPS reconstruction method differs significantly from traditional reconstruction methods, particularly with respect to quality evaluation standards. The similarity between the reconstructed image obtained by the traditional reconstruction approach and the original image is evaluated by the ratio of voxels with the same state value at the same position as the total voxel. A higher ratio indicates better reconstruction similarity. However, this criterion does not apply to MPS reconstruction. The MPS simulation doesn't rely on an exact match between the simulated and training images during reconstruction. The primary goal is to extract details about the training image's microstructure and replicate them in the digital rock that has been reconstructed. Multiple results can be obtained using the MPS simulation, and all of these results can reflect the microstructure characteristics of the training image. The quality of the reconstruction results depends on whether these microstructural features can be copied to the reconstructed digital rock.

Two-point correlation function
When the spatial microstructure changes, the two-point correlation function can show how the phases are related and how they vary and can be used as the evaluation basis of the reconstruction. If two images have the same trend in the twopoint correlation function, this indicates that the attribute phases have similar microstructural features, and vice versa. For the evaluation of multi-mineral shale digital rocks, the two-point correlation function can be defined as the possibility that two voxel positions with a spacing h in the digital rock image demonstrate the same phase state and can be expressed as (Hashemi et al., 2014): where S j is the j phase two-point correlation function, such as the pore, matrix, organic matter, or pyrite phase, f j (u) and f j (u + h) is the phase function of the j phase at position vector u and u + h, respectively, h represents the positional vector.

Pore structure analysis
In the context of statistical analysis, the two-point correlation function just reflects the similarity degree between the training and the reconstruction images.
However, it usually occurs when the low-order statistical functions of each phase in the training image are identical but the morphological microstructure is very different. From a practical application point of view, the reconstructed digital rocks will be adopted to investigate their reservoir physical properties, such as flow features and conductivity. These physical properties are closely correlated to the spatial features, pore structure, distribution, and different phases. Therefore, we need to evaluate the reconstructed digital rocks by analyzing the morphological characteristics of the different attribute phases.
Due to the inherent differences in scale between 2D and 3D images, different morphological feature parameters must be chosen in order to probe the discrepancy between the reconstructed model and the training model used as a benchmark. For comparison and evaluation with 2D images, the percentage of different attribute phases, number, and variety of pore characteristics such as average radius, tortuosity, and fractal dimension were analyzed as the criteria for morphological feature analysis (Maurer et al., 2003;Müllner et al., 2016). Similarly, for comparison and evaluation with 3D images, the percentages of different attribute phases, radius distribution curves (Dong and Blunt, 2009), fractal dimension, tortuosity, and coordination number distribution were used for comparison and evaluation with 3D images.
The percentage of different attribute phases is the proportion of voxels in the overall space of the target attribute phase. Porosity, for example, is the percentage of the pore phase. The radii of the different attribute phases are defined here as the equivalent circle or sphere radius of the Schematic showing the pore radius location of the 2D image, pure white shown in solid pixels, and different depths of blue represent the process of locating the pore radius.

Frontiers in Earth Science
frontiersin.org target attribute phase, as shown in Figure 3, to solve the pore radius of the 2D image. The number of equivalent circles or spheres is the number of target attribute phases, and the distribution curve of the radii of the whole number corresponds to the target attribute phase radius distribution curve. Fractal dimensions can be adopted to show the spatial topology of different attribute phases. The size of the target attribute phase is self-similar to that of the interval with (λ min , λ max ), and the higher and lower limits are the maximum and minimum pore sizes of the target attribute phase, respectively. The relationship between the percentage of the target attribute phase ϕ p and the fractal dimension D f as well as the maximum and minimum dimensions is (Sakhaee-Pour and Xia et al., 2019): where D e is the Euclidean dimension, which is two and three for 2D and 3D space, respectively. The fractal dimension ranges from one to two for 2D space and from two to three for 3D space. The tortuosity, as an important parameter describing the seepage path, depends mainly on material transport properties, including permeability, electrical and thermal conductivity in porous rocks, and is calculated as the proportion of the real length of the seepage channel to the apparent length through the porous medium (Latour et al., 1995;Fu et al., 2021). The coordination number, which is also an important parameter influencing the seepage law, refers to the number of throats connecting the pore to the surrounding pores.

Flow characteristic
Investigating fluid flow in digital rocks, as well as the coexistence and distribution pattern of multi-phase fluids, will shed light on the relationship between a rock's morphological structure and its flow characteristics and provide a solid scientific foundation for studies of relevant applications. The image-based microscopic flow simulation method was developed based on digital rocks combined with computational fluid dynamics methods, which is a common method to simulate the unidirectional and two-phase seepage characteristics of complex rock structures.
For single-phase fluid flow, the finite volume method was utilized to calculate the incompressible Navier-Stokes Equation to determine the velocity field of the digital rocks (Aziz et al., 2018), and the permeability of the digital rocks was calculated. In the flow direction, the inlet and outlet are imposed constant pressure boundary conditions. The fluid-solid interface is a non-slip boundary.
For water-oil two-phase flow, the entire experimental process included oil and water immersion. To simulate oil accumulation in the reservoir, the digital rock was first saturated with water in a drainage process followed by an oil imbibition process. The pore was completely occupied by water after the drainage due to its hydrophilicity. After oil imbibition, the non-wetting oil phase was usually in the pore center, while the wetting water phase preferred to stay at the corners. Then water was injected from the inlet to displace the oil, which aimed to imitate the water flooding process in the realistic reservoir. The simulation approach described above led to the determination of the relative permeability of phases (Soulaine et al., 2021;Yang et al., 2021).

Results and discussion
Four different types of reconstructed simulation cases were adopted to illustrate the practical application capability of the developed algorithm. The two-point correlation function was viewed as a statistical indicator. To evaluate the reconstructed reliability of images and validate the accuracy of the proposed method, the training (reference) image and the reconstructed models were compared with morphological pore structure and fluid flow parameters incorporating absolute and relative permeability curves.

Reconstruction of 2D shale digital rocks from 2D segmented images
The reconstruction of 2D shale digital rocks based on 2D segmented images was chosen as the first case of reconstruction simulation in this study because it is the basis for the subsequent construction of 3D models based on 2D images. In this simulation case, a 2D shale segmentation image of 1000×1000 pixel size was selected as the training image, as illustrated in Figure 4A), which was scanned at a resolution of 60 nm/voxel and generated through combining the proposed algorithm to obtain realizations 1 and 2, as depicted in Figures 4B, C. The images contained organic and inorganic pores, organic matter, matrix, and pyrite. Visually comparing the reconstructed images and the training image, it was found that there was a good similarity between different phases, both in terms of their appearance and shapes. Despite the visual similarity, the visual perception was relatively subjective and could not be quantified for comparison. Therefore, for quantitative comparison, three parameters were selected, namely, the two-point correlation function, morphological structural characteristics, and fluid flow characteristics.
The two-point correlation function was plotted for the evaluation of the statistical change characteristics of the topological relationship between the training image and the reconstructed model, as plotted in Figure 5. The x-axis indicates the hysteresis distance, and the y-axis depicts the  1 and 2, respectively, which were reconstructed using the suggested model. Red indicates the organic pore phase, green shows the organic matter phase, purple indicates the inorganic pore phase, yellow indicates the matrix, and light blue indicates the pyrite phase.

Frontiers in Earth Science
frontiersin.org 08  Frontiers in Earth Science frontiersin.org calculated value of the two-point correlation function of the training image and different reconstructed models of organic and inorganic pores, organic matter phases, and pyrite phases, which show a high degree of agreement, as reflected by the corresponding values of the hysteresis distance and two-point correlation function when the curve ripples and begins to smooth. In a statistical perspective, the two-point correlation function may not be a reliable indicator of spatial morphology and geometric structure, and spatial pore structures with similar two-point correlation functions can also have very different morphological characteristics. For this reason, the relevant morphological structural characteristics are summarized in Table 1. From the percentage of different phases, the values of the percentage of organic pores, organic matter, and inorganic pores in realizations one and two and the training image are very similar. The error ratios are calculated using Eq. 6 to quantify the difference between the reconstructed and training models in terms of the percentage of different phases. Compared with the training images, the error ratios of the percentages of organic pores, organic matter, and inorganic pores in realization one are 1.35%, 3.24%, and 3.05%, respectively, and 4.05%, 3.09%, and 0.7%, respectively, in realization 2. The values of the percentages of matrix and pyrite in realizations one and two are different from those in the training image. Compared with the training image, the error ratios of the percentages of matrix and pyrite in realization one were 6.7% and 40.02%, respectively, and the error ratios of the percentages of matrix and pyrite in realization two were 5.52% and 31.08%, respectively.
where δ i represents the reconstructed model error relative to the training image in the i phase, ϕ ti i is the proportion of the i phase in the training image, and ϕ re i indicates the proportion of the i phase proportion in the reconstructed model.
Considering the number of equivalent regular structures, average radii, and fractal dimensions of organic and inorganic pores, organic matter, matrix, and pyrite, it can be seen that the numerical differences between realizations one and two and the training image are not significant, which indicates that the reconstructed models and the training image exhibit similar morphological structural characteristics. Realizations one and two agree well with the training image in relation to the tortuosity of organic and inorganic pores in x and y-directions, which shows that the reconstructed model and the training image have similar connectivity. The above analysis further confirms the reliability of the proposed algorithm. Because the 2D images involved in this simulation example are mostly discrete, unconnected organic and inorganic pore spaces, the simulation of unidirectional and two-phase fluid flow cannot be carried out, so the permeability and relative permeability cannot be calculated. Therefore, a comparison of the relevant fluid flow characteristics could not be carried out in this simulation example.

Reconstruction of 3D shale digital rocks from 2D segmented images
This study further developed the method of constructing 3D shale digital rocks based on 2D segmented images based on the reconstruction of 2D shale digital rocks based on 2D segmented images to evaluate the algorithm's practical application capability. Figure 6A and B show a real reference model with a pixel size of 400 × 400 × 400 used in this simulation and a cross-Section 2D shale segmentation image selected as the training image, which was scanned at a resolution of 4 nm/voxel and reconstructed by combining the algorithm proposed in this paper to obtain realizations one and two Figure 6C and D. From the figures, it is visible that the images are made of organic pores, organic matter, matrix, and pyrite. Moreover, comparing the reconstructed images with the reference image shows that there is some noise in the reconstructed models, which is unavoidable in the stochastic method. However, there is a lot similarity between the different phases in terms of their appearance and shape. Despite the visual similarity, visual perception is relatively Comparison of (A) the 2D training image (B) 3D reference image and (C) and (D) realizations one and two which are reconstructed from the training image using the proposed approach. Black indicates the organic pores, light blue indicates the organic matter, yellow indicates the matrix, and pink represents the pyrite.

Frontiers in Earth Science
frontiersin.org subjective and cannot be compared quantitatively; thus, twopoint correlation function, morphological structure features, and fluid flow features are selected for quantitative comparison. Figure 7 shows the change in the two-point correlation function with distance for different phases. The results demonstrate a high degree of consistency, which is reflected in the corresponding values of the hysteresis distance and two-point correlation function when the curve ripples and begins to smooth. When the curve is smoothed, the two-point correlation functions of the matrix have different values compared to the ones of organic pores and organic matter. This discrepancy can be traced back to a proportional change in the matrix between the reference image and the reconstructed versions of the image. However, the consistent trends of both curves indicate that the statistical characteristics of the matrix are similar.
Spatial pore structures with similar two-point correlation functions may have very different morphological characteristics. Table 2 calculates the petrophysical properties and pore structure parameters. The percentage values of the organic pores and organic matter in realizations 1 and 2, and the reference image were very close. Furthermore, the percentage values of matrix and pyrite in realizations 1 and 2, and reference image were somewhat different. In order to further quantify the percentage values of the various phase differences in the considered images, the error ratios were calculated using Eq. 6. Compared with the reference image, the error ratios of the occupancy ratios of organic pores and organic matter in realization one were 11.10% and 6.80%, respectively, and 9.13% and 3.25% in realization 2. On the other hand, the error ratios of the percentages of matrix and pyrite in realization one were 20.65% and 39.66%, respectively, and 23.49% and 54.86% in realization 2. Despite the significant differences in the ratios of matrix and pyrite in the reconstructed models, the pore space parameters of organic pores and organic matter have the most impact on the flow.
Comparing the equivalent rule radius distributions of organic pores, organic matter, matrix, and pyrite, as depicted in Figure 8, the distribution curves of both realizations one and two are close to those of the reference image. This demonstrates that the reference and reconstructed images have similarity in morphological structure characteristics, which can also be confirmed by the comparison of the fractal dimensions in Table 2. The distribution curves of the coordination numbers of the models are calculated and plotted, as shown in Figure 8B, showing a similar distribution, which indicates that the reference and reconstructed images have similar connectivity. The above analysis provides proof that the proposed algorithm is reliable.

FIGURE 7
Comparison of two-point correlation functions of reference image and reconstructed images for different phases(A) organic pore phase (B) organic phase, (C) matrix phase Black represents the training image, red represents realization one and blue represents realization 2.

Frontiers in Earth Science frontiersin.org
An examination of the two-point correlation function and morphological structure features shows that the proposed algorithms are accurate. However, both parameters are based on small amounts statistical data, and the prediction of absolute and relative permeability based on fluid flow simulation laws is the key to assessing the quality of the reconstructed images. Table 2 shows that the difference between the absolute permeability of the reconstructed images and that of the reference image is small, which further verifies the reliability of the reconstructed model. Compared with the testing of absolute permeability, the testing of the relative permeability is quite difficult because it is highly sensitive to the morphological structure and 3D connectivity. By plotting the relative permeability of the reconstructed models and the reference model, one can note that for the relative permeability of water, the reference and reconstructed model curves follow the same trend and almost overlap (see Figure 9). However, for oil-phase permeability, they have some differences, but their basic curve trends are similar. Comparing the values of connate water saturation, residual oil saturation, and iso-permeability point, the reconstructed models were close to the reference model.

Reconstruction of 2D shale digital rocks from 2D gray images
The simulation cases illustrate how the proposed algorithm is used in this study for multi-mineral digital rock reconstruction with segmented images. In the

FIGURE 8
Comparison of pore structure distribution of the reference image and two realizations for different phases. (A) Distribution of organic pore radius, (B) distribution of organic pore coordination, (C) distribution of organic matter radius, (D) distribution of matrix radius, and (E) distribution of pyrite radius.

Frontiers in Earth Science
frontiersin.org simulated instances below, we further exhibit the capability of the algorithm in continuous-phase porous media reconstruction applications. As shown in Figure 10A, a 2D shale gray image of 1000×1000 pixel size was chosen as the training image. It was scanned at a resolution of 60 nm/voxel, combined with the proposed algorithm for obtaining realizations one and 2 (see Figure 10B and C). Figure 11 shows the gray value distribution curves of the reconstructed and the training images. It can be shown that the gray value distribution curves of the reconstructed and the training images show a similar three-peak distribution, especially the performance at the peak, which indicates that the two have similar texture characteristics. Based on the delineated gray value as the threshold, the gray image is segmented into different organic pores, organic matter, inorganic matter pores, matrix, and pyrite (see Figure 12). A visual comparison between the reconstructed models and the training image revealed a good similarity between the different phases, both in terms of their appearance and shape. Although the visual perception is relatively subjective and cannot be quantified, two-point correlation functions, morphological structure features, and fluid flow characteristics were selected for quantitative comparison.
To evaluate the variation features of the spatial distribution on the reconstructed and training images for different phases, the change of the two-point correlation function was plotted ( Figure 13). As evidenced by the corresponding values of the lag distance when the curves ripple trend and the points start to smooth out, the two-point correlation functions of the training image and various reconstructed models of different phases such as the organic pore, organic matter, inorganic pore, and pyrite phases exhibit a high degree of agreement.
The reconstructed models were eventually used to simulate the fluid flow process and predict the seepage properties, which mainly depend on the morphological structure characteristics computed in Table 3. In terms of the percentage of different phases, the percentages of different phases in realizations one and two and the training image are very close to each other. Compared with the training image, the error ratios of the percentages of the organic pore, organic matter, and inorganic pore phases in realization one were 0.68%, 4.56%, and 4.23%, respectively, and 3.38%, 3.82%, and 0.59%, respectively, in realization 2. The values of the percentages of matrix and pyrite in realizations one and two and the training image are somewhat different. Compared with the training image, the error ratios of the percentages of matrix and pyrite in realization one are 6.25% and 36.84%, respectively, and the error ratios of the percentages of matrix and pyrite in realization two are 3.96% and 20.55%, respectively.
Considering the number of equivalent regular structures, mean radii, and fractal dimensions of phases, it can be understood that the numerical differences between realizations one and two and the training image are not

FIGURE 9
Comparison of the predicted relative permeability during the imbibition displacement. Black represents the computed relative permeability of reference, and red and blue are the computed relative permeabilities of realizations 1 and 2, respectively. The solid line is the water relative permeability, and the dashed line is the oil relative permeability.

FIGURE 10
Comparison of 2D shale gray digital rocks (A) Training image, (B) realization 1and (C) realization 2. Realizations are reconstructed based on the training image using the proposed approach.

Frontiers in Earth Science
frontiersin.org significant, which indicates that the reconstructed model and the training image have similar morphological structure characteristics. The tortuosity of Realizations one and two match well with the training image for organic and inorganic pores in x and y directions, which shows that the reconstructed model and the training images have similar connectivity. The above analysis further confirms the reliability of the proposed algorithm. Unidirectional and two-phase fluid flow simulations cannot be performed in this simulation example because the majority of the 2D images are discontinuous, disconnected organic and inorganic pore spaces. As a result, the permeability and relative permeability cannot be obtained. Therefore, a comparison of the relevant fluid flow characteristics is not carried out in this simulation example.

Reconstruction of 3D shale digital rocks from 2D gray images
The study further develops the construction of 3D models based on 2D gray-scale images to evaluate the practical application capability of the proposed algorithm. In this simulation case, a 400×400×400 pixel size is chosen as the real reference model, as shown in Figure 14B, from which a cross-sectional 2D shale gray image is chosen as the training image (see Figure 14A), which was scanned at a resolution of 4 nm/voxel, and reconstructed by combining the proposed algorithm in this paper to obtain realizations 1 and 2, as depicted in Figure 14C and D. Figure 15 shows the gray value distribution curves of the reconstructed and the reference images. It can be shown that the

Contents
Training image Realization#1 Realization#2 The petrophysical properties Organic pores proportion/% gray value distribution curves of the reconstructed model image and the reference image show a similar multi-peak distribution, especially the performance at the peak, which indicates that the two have similar texture characteristics. Based on the demarcated gray values as the thresholds, the gray image is segmented into different organic and inorganic pores, organic matter, matrix, and pyrite (see Figure 16). From the figure, it can be shown that the image contains organic pores, organic matter, matrix, and pyrite and visually compares the reconstructed models with the reference image. However, there is considerable noise in the reconstructed model, which is inevitable with the stochastic method. There is a good resemblance between the various phases in terms of their appearance and shape. However, it cannot be quantitatively compared since it is so subjective. As a result, three types of indicators were chosen for quantitative comparison: two-point correlation functions, morphological pore structure features, and fluid flow features.
Here again, the two-point correlation functions of the reference images of organic pores (see Figure 17A) and organic matter (see Figure 17B) of different reconstructed models and the reference image show a high degree of consistency, which is reflected when the curve ripples and begins to smooth. Compared to the organic pore and organic matter, the two-point correlation function of the matrix (see Figure 17C) has different values when it is smoothed, which is related to the difference in pyrite in the reconstructed images. This is related to the difference in the proportion of pyrite in the reconstructed image and the reference model. However, the consistent trends of the two curves indicate that the statistical characteristics of pyrite are similar.
The relevant pore structure parameters are computed in Table 4, and the occupancy ratio values of the organic pore and organic matter in realizations one and two and the reference image are very close. For further quantitative measurements of the different phase occupancy ratios in the reconstructed and the reference models, differences were evaluated by calculating the error ratios defined in Eq. 6. Compared with the reference image, the error ratios of organic pores and organic matter occupancy ratios in realization one were 2.13% and 0.82%, respectively, and 0.008% and 2.62%, respectively, in realization 2. The values of the percentages of matrix and pyrite in realizations one and two and the reference image were somewhat different. Compared with the reference image, the error ratios of the percentage of matrix and pyrite in realization one were 12.37% and 35.84%, respectively, and the error ratios of the percentages of matrix and pyrite in reconstructed realization two were 18.72% and 49.56%, respectively. In spite of the significant disparities in the reconstructed models' percentages of pyrite and matrix, the pore structure characteristics had the biggest impact on the fluid flow characteristics.

FIGURE 11
Comparison of the gray value distribution of the training image and different realizations based on the proposed method. Black is the gray value distribution of the training image. Red and blue are the gray value distributions of realizations 1 and 2, respectively. The green dotted line is the gray value used to segment 2D shale gray digital rocks to obtain 2D shale segmented digital rocks.

FIGURE 12
2D shale segmented digital rocks (A) training image, (B) and (C) are realizations 1 and 2, respectively, which are obtained by segmenting the reconstructed gray images. Red indicates the organic pore phase; green indicates the organic matter phase, purple indicates the inorganic pore phase, yellow indicates the matrix, and light blue indicates the pyrite phase.

Frontiers in Earth Science frontiersin.org
Comparing the equivalent rule radius distributions of organic pores, organic matter, matrix, and pyrite, as depicted in Figure 18, the distribution curves of both realizations one and two are close to those of the reference image. This means that the reconstructed images and the reference image show similar morphological and structural features, which can also be confirmed by the comparison of fractal dimensions in Table 4. Furthermore, the plot of their coordination numbers reveals similarity in distribution, which means that they have similar connectivity (see Figure 18). The above analysis further confirms the reliability of the proposed algorithm.
Once again, the assessment of the two-point correlation function and morphological structure features shows that the proposed algorithms are accurate. However, as previously stated, both have the disadvantage of being limited to small data sets. Thus, the prediction of absolute and relative permeability based on fluid flow simulation laws is the key to assessing the quality of the reconstructed images. The results in Table 4 reveal that the difference between the absolute permeability of the reconstructed images and that of the reference image is slight, further validating the accuracy of the proposed approach. Testing the oil-water relative permeability is difficult, because it is highly sensitive to the morphological structure and 3D connectivity. By plotting the relative permeability of the reconstructed images and the reference image, as depicted in Figure 19, it can be shown that for the relative permeability of water, their curves follow a similar trend and almost overlap, but for oil relative permeability, the reconstructed models differ slightly from the reference model. However, their basic trends are similar. Compared to the values of the connate water saturation, residual oil saturation, and iso-permeability point, the reconstructed models are close to the reference model.

Conclusion
Given that shale oil and gas reservoirs develop multiple types of complex pores and minerals and exhibit multi-scale characteristics, the most widely used experimental physical methods are limited by the contradiction of resolution and representative elementary volume, which makes it difficult to accurately construct shale multi-mineral digital rocks and thus carry out subsequent fluid flow simulation studies. Therefore, it is of utmost importance to develop a numerical approach to construct shale multi-mineral digital rocks.

FIGURE 13
Different phases two-point correlation functions (A) Organic pore, (B) organic, (C) inorganic pore phase, and (D) pyrite phase. It should be noted that the images were obtained by segmenting the reconstructed gray images.

Frontiers in Earth Science
frontiersin.org 1 This paper proposes an innovative approach for constructing shale multi-mineral digital rocks from 2D images (segmented or gray), which is based on multi-point statistics, and can realize the construction of segmented and gray phase shale multi-mineral digital rocks in 2D and 3D. This work shows that it is possible to obtain 3D shale multi-mineral digital rocks using only one 2D training image as the input, combined with the proposed algorithm, while reproducing the pore structures in a real model. 2 In different simulation cases, the two-point correlation function and morphological pore space features of each

FIGURE 14
Comparison of 3D shale digital rocks (A) Training image, (B) Reference model, (C) and (D) are realizations one and two which are reconstructed based on the training image using the proposed method. Red represents the organic pore phase; green represents the organic phase, purple represents the inorganic pore phase, yellow represents the matrix phase, and light blue represents the pyrite phase.

FIGURE 15
Comparison of the gray value distribution of the reference image and different realizations based on the proposed approach. The green dotted line is the gray value used to segment 3D shale gray digital rocks to obtain 3D shale segmented digital rocks.
Frontiers in Earth Science frontiersin.org

FIGURE 16
3D shale segmented digital rocks (A) Reference image, (B) and (C) are realizations one and two which are obtained by segmenting the reconstructed gray images. Black represents the organic pore phase, light blue represents the organic phase, yellow represents the matrix phase, and pink represents the pyrite phase.

FIGURE 17
Two-point correlation functions of different phases for the training image and two realizations obtained by segmenting the reconstructed gray images, (A), (B), and (C) are comparisons of two-point correlation functions of the reference image and different realizations (realization one and 2) obtained by segmenting the reconstructed gray images for the organic pore, organic matter, and matrix phase, respectively. Black represents the training image, red represents realization one and blue represents realization 2.
Frontiers in Earth Science frontiersin.org 18 phase are well reproduced in the process of constructing 2D shale digital rocks based on 2D training images, both in the segmented and gray regions. However, for the process of constructing 3D shale digital rocks based on 2D images, the two-point correlation function and morphological structure of organic pores and organic matter are well reproduced, both in the segmented and gray, unfortunately, the matrix and pyrite are relatively poorly reconstructed. Mathematically, the construction of 3D images from 2D images can be viewed as a spatially inverse problem that only borrows the pore structure from 2D images and restores the pore structure in 3D images, which has a large amount of randomness and uncertainty. However, based on the performance of absolute and relative permeability, the values and curve distributions of the reconstructed images are similar to those of the reference image, which indicates that the reconstructed model reproduces the microstructure of the real 3D space. Improving the accuracy of the algorithm is a critical task for moving the research forward in the next step.

FIGURE 18
Comparison of pore structure distribution of the reference model and two realizations reconstructed by the proposed approach. (A) Organic pore radius distribution, (B) organic pore coordination distribution, (C) organic radius distribution, (D) matrix radius, and (E) pyrite radius. In the pictures (A), (C), (D), and (E), black is the training image distribution, and red and blue are the distributions of realizations 1 and 2, respectively. In (B), green is the training image distribution, and red and blue are the distributions of realizations 1 and 2, respectively.

Frontiers in Earth Science
frontiersin.org 3 The algorithm proposed in this paper provides a straightforward and fast tool for 3D multi-mineral reconstruction of 2D slices. The application of the proposed algorithm through multiple simulation cases shows that the algorithm is reliable and time-and costeffective compared to experimental methods. Although the algorithm is proposed in the context of shale, but the algorithm can be extended to other types of reservoirs as well.

Data availability statement
The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.

Author contributions
LL, HS, and JY conceived the project; LL coded the related software and carried out the simulation and analyses and wrote the original manuscript; GI suggested the organization of the manuscript; HS and LZ provided technical advice throughout the project; YY and KZ assisted with the writing of this paper; JY supervised the overall project, and all author contributed to revising and improving the manuscript.

FIGURE 19
Comparison of predicted relative permeabilities to the oil and water phases during imbibition process. Black represents the computed relative permeability of the reference model, red and blue are the computed relative permeabilities of realizations 1 and 2, respectively. The solid line is the relative permeability of the water phase, and the dashed line is the relative permeability of the oil phase.
Frontiers in Earth Science frontiersin.org