# Challenges and Prospects of Digital Core-Reconstruction Research

Authors: Linqi Zhu; Chong Zhang; Chaomo Zhang; Xueqing Zhou; Zhansong Zhang; Xin Nie; Weinan Liu; Boyuan Zhu
Journal: Geofluids
 (2019)

---

## Abstract

The simulation of various rock properties based on three-dimensional digital cores plays an increasingly important role in oil and gas exploration and development. The accuracy of 3D digital core reconstruction is important for determining rock properties. In this paper, existing 3D digital core-reconstruction methods are divided into two categories: 3D digital cores based on physical experiments and 3D digital core stochastic reconstructions based on two-dimensional (2D) slices. Additionally, 2D slice-based digital core stochastic reconstruction techniques are classified into four types: a stochastic reconstruction method based on 2D slice mathematical-feature statistical constraints, a stochastic reconstruction method based on statistical constraints that are related to 2D slice morphological characteristics, a physics process-based stochastic reconstruction method, and a hybrid stochastic reconstruction method. The progress related to these various stochastic reconstruction methods, the characteristics of constructed 3D digital cores, and the potential of these methods are analysed and discussed in detail. Finally, reasonable prospects are presented based on the current state of this research area. Currently, studies on digital core reconstruction, especially for the 3D digital core stochastic reconstruction method based on 2D slices, are still very rough, and much room for improvement remains. In particular, we emphasize the importance of evaluating functions, multiscale 3D digital cores, multicomponent 3D digital cores, and disciplinary intersection methods in the 3D construction of digital cores. These four directions should provide focus, alongside challenges, for this research area in the future. This review provides important insights into 3D digital core reconstruction.


---

## Body

## 1. Introduction

With continuous exploration and development in the petroleum industry, unconventional reservoirs, such as tight oil and gas, shale oil and gas, coalbed methane, and natural gas-hydrate reservoirs, have received more attention than conventional oil and gas reservoirs; thus, these unconventional reservoirs are becoming key areas for exploration and development [1–10]. Unconventional reservoirs are quite different from conventional reservoirs; these differences are reflected mainly in the microgeological characteristics of rocks [11, 12]. The rocks in unconventional reservoirs are characterized by various pore types, complex pore structures, diverse mineral compositions, and complicated oil and gas distributions [13–15]. Complex geological features result in complex seepage, electrical, acoustic, elastic, radioactive, and nuclear magnetic resonance (NMR) features of the rock [16–19]. Consequently, exploration methods that are widely used for conventional reservoirs are difficult to apply directly to unconventional reservoirs, which require targeted research [20–25]. Conducting an effective study on unconventional reservoirs is difficult because of the lack of accurate data, the high cost, and the need to improve the core experimental methods that are employed for the characterization of unconventional reservoir rocks. Since the concept of a digital rock-physics experiment was proposed, research on virtual rock physics based on digital cores has constituted a popular research topic coincident with the development of related disciplines and great improvements in computational capabilities [26–32]. Virtual rocks boast a number of advantages; for example, these virtual rocks are intuitive, they are strongly reproducible and nondestructive in experiments, they have easily adjustable parameters, and they are relatively inexpensive and require only short experimental periods [33]. Therefore, virtual rock physics has great potential as an alternative experimental method, especially for shale reservoirs. Virtual rock physics also represents an important development direction in petrophysics.For virtual rock physics research that uses digital cores, guaranteeing the accuracy in a simulation experiment is based on the accuracy of modelling the numerical core. The oldest digital core models are the abstract capillary bundle model and the random pore network model; however, the concept of digital cores was unknown when these models were initially created. The capillary bundle model approximates the pore space of the rock by using capillary tubes of varying radii [34–36]. Although this method is extremely abstract with respect to the pore network of a rock, it is valuable for evaluating the seepage characteristics of conventional reservoirs and primary tight-pored sandstone reservoirs [37, 38]. In contrast, the random pore network model uses pores and pore throats to form a mesh structure from real experimental methods, such as mercury intrusion experiments or NMR experiments, to simulate rocks with less complex structures [39]. Nevertheless, the abovementioned numerical core models simplify the pore space of the rock, so accurately reflecting the characteristics of unconventional reservoirs is difficult.The random pore structure model simplifies only the pore space of the reservoir rock and thus cannot fully reveal the pore structure of a complex reservoir. For example, characterizing the full pore size distributions of shale reservoirs through pore structure experiments is difficult because overly complex pore structures can reduce the experimental porosity measurement accuracy. In addition, the current pore structure model cannot reflect the influences of the rock matrix on the physical properties of the reservoir; hence, only the fluid-transfer characteristics can be studied. Greyscale images of rocks can be obtained with the application of X-ray computed tomography (CT) scanning and other techniques to the study of reservoir rock microstructures. The different greyscale values in these images represent the different components of the rock. An algorithm can be used to divide the obtained three-dimensional (3D) gradation into two components, namely, the rock matrix and the rock pores, and different integer values are used to represent different components. In this way, a digital rock, that is, a modern digital core, is obtained. Compared to the previous pore network model, a digital core can both reflect the complex pore space of a reservoir rock and display its solid skeleton [40, 41].At present, the digital cores that are used to study unconventional reservoirs are usually derived from rock-based digital images. Numerical simulations of physical properties based on digital cores are called digital petrophysical experiments; among them, digital rock physics experiments of acoustic, electrical, NMR, and seepage properties are currently available. The real structure of the rock is obtained in this fashion, so the numerical simulation results of the physical properties of the corresponding rock match the experimental measurements.Since the inception of this new field, many scholars have studied the reconstruction of 3D digital cores [42]. According to the above description, digital cores are a very important component of modern petrophysical research. At present, 3D digital cores can be reconstructed in many ways, that is, by various experimental methods and algorithms, and the accuracy of the 3D digital core directly affects the implications of subsequent simulation experiments. Therefore, we must review the progress of 3D digital core reconstruction and provide a reasonable outlook for its efficacy in future studies.


## 2. Status of Utilizing Physical Experimental Methods to Reconstruct 3D Digital Cores

Physical methods are employed to directly reconstruct 3D digital cores; that is, the real structure of the rock samples is obtained by 3D scanning or continuous slice scanning with instruments and various other physical means. Three main physical experimental methods are used to create 3D digital cores: confocal laser scanning, serial section imaging, and X-ray CT scanning. Sections2.1 and 2.2 of this paper introduce the principles, resolution, advantages, and disadvantages of confocal laser scanning, serial section imaging, and X-ray CT scanning.


### 2.1. Confocal Laser Scanning

Confocal laser scanning is used to obtain the 3D pore distribution of a sample by using a laser scanning confocal microscope (LSCM) [43, 44]. The sample must be prepared appropriately before performing laser scanning, and the thickness of the sample must be controlled during the preparation of the sample. This control is required because the limited maximum penetration depth of the LSCM is approximately 100 μm, which makes this method unsuitable for thicker samples. After drying, the sample is injected into dyed epoxy by using a vacuum and pressurized infusion. The epoxy resin fluoresces when excited by the laser, and the LSCM detects this fluorescence. Thus, the LSCM can describe the 3D pore distribution of a rock sample.Laser scanning confocal microscopy [45] uses a laser beam to illuminate the sample; the laser beam is then reflected by the beamsplitter and focused on the sample, and each point in the sample is scanned. If an excited dye is present in the sample, then the fluorescence is returned to the beamsplitter through the incident light path, after which it is focused and collected by a photomultiplier tube (PMT) and displayed on a computer (in greyscale). The detection area of the microscope can be precisely controlled by the computer, moved along the surface of the sample, and probed into different depth regions of the sample to obtain a digital core. Generally, the resolution of confocal laser scanning can reach the micrometre-submicron scale.However, this method has many drawbacks. Foremost among them, the obtained digital core cannot be overly thick; that is, the sample usually cannot exceed a thickness of 200μm (Figure 1). In addition, the penetration depth is generally less than 1 mm. The second weakness of this technique is that it requires the injection of an epoxy resin that contains a fluorescent agent. This requirement constitutes a substantial problem for three reasons. First, the injection of epoxy resin directly destroys the core, so the core can no longer be used for other purposes. Second, epoxy cannot be injected into isolated voids; thus, ineffective pores within the rock are undetectable. Third, as a liquid, epoxy resin has a certain viscosity and is difficult to inject into excessively small pores, such as organic pores. These drawbacks can prevent the identification of micropores regardless of the accuracy of the probing instrument. The above analysis shows that the potential for further development of the LSCM method is inadequate, and this approach may be replaced by other methods in the future. Accordingly, confocal laser scanning is not widely used in practical research.Figure 1
Reconstruction of 3D digital cores using confocal laser scanning method (the pore space is opaque (colored) and the solid matrix is translucent). It is obvious that the thickness of 3D digital cores obtained by LSCM method is low. This leads to too little information on the obtained 3D digital core, and the results of numerical simulation using 3D digital cores are not representative. Figure from [43].


### 2.2. Imaging Methods Based on a Serial Section Rock Overlay

#### 2.2.1. Ordinary Sequence Two-Dimensional (2D) Slice Superposition Method

Rock imaging methods based on a serial section rock overlay are destructive imaging techniques, of which two exist: the ordinary sequence 2D slice superposition method and the focused ion beam scanning electron microscopy (FIB-SEM) physical imaging method [46, 47]. In the former, the rock sample is prepared by polishing it to obtain a relatively flat plane, after which a high-magnification microscope is used to photograph the polished rock sample surface to obtain a microscopic image of the core surface. Subsequently, a layer of rock specimen is cut parallel to the polishing surface, and the cut rock specimen is further polished and photographed with a high-magnification microscope. This process is repeated until a 3D digital core with the desired thickness is obtained. Finally, experimental images of the core are obtained by combining these photographs. This method yields nanoscale and high-resolution core images [48–50].


#### 2.2.2. FIB-SEM Physical Imaging Method

The ordinary sequence 2D slice superposition method can obtain a higher-resolution 3D digital core, but this technique uses an electron beam to polish the sample surface, which produces static electricity at the surface; therefore, imaging the surface is not beneficial. As an alternative, FIB-SEM was first developed in 1988. The FIB-SEM system utilizes two beams and is thus referred to as the dual-beam system [51]. FIB-SEM can be simply understood as a single-beam FIB combined with an SEM. Differently from the ordinary-sequence 2D slice superposition method, the use of ion-beam polishing has the advantage of producing less static electricity and better imaging quality. Therefore, the FIB-SEM imaging method is more applicable to this study than the ordinary sequence 2D slice superposition method.Although the resolution can reach the nanoscale based on the rock serial section overlay imaging method, the FIB-SEM technique is problematic because it destroys rock samples and has a slow modelling speed (Figure2). For a sample with a surface area of 50μm×50μm, the FIB abrasion of a sample with a thickness of 0.1 μm can take minutes, and only 5-20 images can be scanned per hour, which substantially slows the modelling process. In addition, the faster the sample erodes, the poorer the polishing effect becomes; thus, the quality of the image is affected. Therefore, increasing the speed of the serial section rock overlay imaging method is difficult [52, 53]. Another disadvantage of this method is the cutting and subsequent destruction of the rock sample. In addition, when a sample is recut, the thickness between adjacent slices is directly reduced during the cutting, and the pore structure of the rock sample is destroyed. Consequently, an inaccurate pore structure is obtained for the rock sample. This phenomenon results in a relatively poor pore structure after splicing that does not completely reflect the original pore structure. These three disadvantages restrict the development of this method. However, the accuracy of this technique is still better than that of focused scanning, and the nanoscale resolution is favourable for the identification of organic pores in shale rocks. Therefore, this method is used more often than focused scanning [54–57].Figure 2
The multicomponent 3D digital core reconstructed by the FIB-SEM method. Different components are reflected in different greyscales. The 3D digital core obtained by this method has a high resolution, but the disadvantage is that the scanned sample will be destroyed.


#### 2.2.3. X-Ray CT Scanning Method

The basic principle of CT, which was proposed by G.N. Hounsfield in 1969 [58], involves the use of X-ray scanning to obtain the thickness of an object as the detector receives the X-rays through the sample. The X-rays are converted into electric signals that consist of visible light through a photoelectric converter and are then converted into digital signals through an analogue-digital converter; finally, the signals are processed by a computer to obtain 3D greyscale volume data. During the actual scanning process, the sample is subdivided into voxels of the same volume, and the X-ray attenuation coefficient of each voxel is determined by scanning and arrangement in a matrix, which is then converted into greyscale volume data. The attenuation coefficient is the key to distinguishing the material components of a digital core. Typically, the number of X-rays that are absorbed by a sample depends on the density of the components in the sample. Therefore, the absorption coefficient reflects the material composition of the sample.The attenuation coefficient of X-rays that pass through a component can be determined by Beer’s law as follows:(1)I=I0⋅exp−μx,where I0 is the incident X-ray intensity, I is the X-ray intensity after attenuation, μ is the linear attenuation coefficient, and x is the length of the X-ray penetration path.If the sample consists of many components (such as rocks with more complex minerals), the above equation becomes the following:(2)I=I0⋅exp∑i−μixi,where i represents the i−th component.In the early 1980s, the first microcomputed tomography system (micro-CT) was developed to increase the vertical resolution to the level of the horizontal resolution [59]. In the early 1990s, Dunsmuir proposed the application of the micro-CT method in the field of petroleum exploration [60]. Later, Coenen et al. used micro-CT techniques to construct submicron digital cores [61]. Currently, two micro-CT systems are commonly employed: a benchtop micro-CT scanning system that generates X-rays with an industrial X-ray generator and a synchronous acceleration micro-CT system that uses a synchrotron as an X-ray generator [62]. Generally, the digital core resolution that is obtained by a synchronous acceleration micro-CT scanning system is higher than that from a desktop micro-CT scanning system, but the resolution of a micro-CT scanning system is typically still within the micrometre to submicrometre range. Figure 3 shows a 3D sandstone core that was constructed from CT scanning by Arns et al. at the Australian National University, and Figure 4 shows the Fontainebleau sandstone that was constructed with CT scanning at Brookhaven National Laboratory.Figure 3
The use of X-ray CT scanning method is the establishment of a 3D digital core. The transparent part is the skeleton; the other part is the pores.Figure 4
The Fontainebleau sandstone dataset with a binary segmentation. The red part is the skeleton; the black part is the pores. The total size of the dataset is288×288×300 and a voxel edge length is 7.5 μm. Figure from [63].In theory, the X-ray resolution limit should be 0.005 nm, but the original micro-CT technique cannot probe at a resolution greater than the submicron scale. The key restriction is the size of the X-ray source and the pixel size of the detector. The detector size of the X-ray source cannot reach the nanometre scale. Fortunately, the resolution of a current radiation source can reach the nanometre scale, and the size of the radiation source can be constrained to achieve the resolution of nanometre-scale CT (nano-CT) [64, 65]. After the appearance of the nano-CT technique, nanometre CT systems based on visible light optical systems, synchrotron radiation sources, and X-ray optical systems have been proposed. For example, Lawrence Berkeley National Laboratory obtained 15 nm resolution CT images. With continuous improvements in CT scanning technology, the resolution of this method can be further improved.Compared to the above two methods, the greatest advantage of the X-ray CT scanning method is that the sample is not damaged during the construction of the digital core. In addition, current nano-CT systems can achieve nanoscale resolution with the same accuracy as FIB-SEM systems. Table1 shows the sample size, instrument resolution, and main application areas of the CT and FIB-SEM analysis techniques.Table 1
Instrument parameters and applications for 3D digital core reconstruction.
Instrument type
Sample size
Instrument resolution
Main application areas
Micro CT
Core column scanning(diameter 15-25 mm)
12-30μm
Geological description, heterogeneity analysis
Core column scanning(diameter 1-25 mm)
0.5-12μm
Fine description, pore network generation
Nano-CT
Less than 1 mm
65-150 nm
For the scanning of complex structural samples such as shale, microcrystal, tight gas rock, and analysis of pore network
FIB-SEM
Less than 1 mm
2-150 nmAmong the abovementioned methods, the X-ray CT scanning technique is currently the most popular. However, this approach also suffers from some problems. For instance, the sample size of a CT experiment is too small. For rocks with a high degree of heterogeneity, the results from CT scanning are not representative, so one must reconstruct the multiscale digital core. Nevertheless, the applicability of CT experiments is generally broader than that of FIB-SEM experiments.


#### 2.2.4. Magnetic Resonance Imaging (MRI) Method

The magnetic resonance imaging (MRI) method was first proposed in 1973 by Lauterbur [66]. Subsequently, MRI methods were gradually applied in the medical and biological communities, among others. In 1995, Xiao applied the MRI method to the field of petroleum engineering for the first time [67], after which this technique was widely employed in both petrophysics and petroleum engineering [68–73]. The use of the MRI method to image fluids in rocks to reflect their structures has constituted an important guide in the visual determination of microscopic pore structures.A nucleus with a nonzero magnetic moment undergoes a phenomenon known as NMR; that is, Zeeman splitting occurs in the spin-energy level of the external magnetic field, and the resonance absorbs radiation that is associated with radio waves at a certain frequency. NMR signals are generated during NMR within the nucleus. These NMR signals are collected by a special device, and a 2D inverse Fourier transform is used to obtain a reconstructed image of a spatial point. The entire target object can be imaged by collecting the signals at each point in the range of the target object.Compared to the LSCM, FIB-SEM, and X-ray CT methods, the MRI method is limited by a very obvious disadvantage: an excessively low imaging resolution [74, 75]. For tight sandstones and shale rocks, obtaining a large amount of microscopic information with the MRI method is difficult. In addition, NMR experiments require that the pores of the rock be saturated with fluids, which is inherently difficult to achieve for very dense rocks. Moreover, the type of fluid within the rock can affect the imaging results [76, 77]. The NMR inversion algorithm is another factor that influences the imaging quality; accordingly, research on NMR inversion algorithms is ongoing [78–80]. In summary, MRI is not currently employed to construct digital cores because of its limitations.Many relatively simple pore network models are also available [81, 82], and experiments such as mercury intrusion and one-dimensional NMR can be used to image the porosity of a rock [83, 84]. However, the former cannot accurately reflect the heterogeneity of rocks. Moreover, the pore shape and overall pore network in a simple pore network model are fixed, and the pores are arranged in a given orientation; these principles are in contrast to the idea that the digital core must reflect the rock features, including the heterogeneity, to the greatest extent possible. Unfortunately, experiments such as mercury intrusion reflect the pore structure only from the side of the sample, which is not intuitive and has multiple solutions. Therefore, the above two methods cannot be called digital core construction methods.Overall, the most promising and suitable physical experimental method that is currently available for the reconstruction of a digital core is CT scanning. The FIB-SEM reconstruction method can also obtain digital cores with relatively high resolution, but the samples will be destroyed in the process. Nevertheless, as the resolution of these instruments continues to increase, the accuracy of physical experimental techniques in the reconstruction of 3D digital cores will also continue to improve.


## 2.1. Confocal Laser Scanning

Confocal laser scanning is used to obtain the 3D pore distribution of a sample by using a laser scanning confocal microscope (LSCM) [43, 44]. The sample must be prepared appropriately before performing laser scanning, and the thickness of the sample must be controlled during the preparation of the sample. This control is required because the limited maximum penetration depth of the LSCM is approximately 100 μm, which makes this method unsuitable for thicker samples. After drying, the sample is injected into dyed epoxy by using a vacuum and pressurized infusion. The epoxy resin fluoresces when excited by the laser, and the LSCM detects this fluorescence. Thus, the LSCM can describe the 3D pore distribution of a rock sample.Laser scanning confocal microscopy [45] uses a laser beam to illuminate the sample; the laser beam is then reflected by the beamsplitter and focused on the sample, and each point in the sample is scanned. If an excited dye is present in the sample, then the fluorescence is returned to the beamsplitter through the incident light path, after which it is focused and collected by a photomultiplier tube (PMT) and displayed on a computer (in greyscale). The detection area of the microscope can be precisely controlled by the computer, moved along the surface of the sample, and probed into different depth regions of the sample to obtain a digital core. Generally, the resolution of confocal laser scanning can reach the micrometre-submicron scale.However, this method has many drawbacks. Foremost among them, the obtained digital core cannot be overly thick; that is, the sample usually cannot exceed a thickness of 200μm (Figure 1). In addition, the penetration depth is generally less than 1 mm. The second weakness of this technique is that it requires the injection of an epoxy resin that contains a fluorescent agent. This requirement constitutes a substantial problem for three reasons. First, the injection of epoxy resin directly destroys the core, so the core can no longer be used for other purposes. Second, epoxy cannot be injected into isolated voids; thus, ineffective pores within the rock are undetectable. Third, as a liquid, epoxy resin has a certain viscosity and is difficult to inject into excessively small pores, such as organic pores. These drawbacks can prevent the identification of micropores regardless of the accuracy of the probing instrument. The above analysis shows that the potential for further development of the LSCM method is inadequate, and this approach may be replaced by other methods in the future. Accordingly, confocal laser scanning is not widely used in practical research.Figure 1
Reconstruction of 3D digital cores using confocal laser scanning method (the pore space is opaque (colored) and the solid matrix is translucent). It is obvious that the thickness of 3D digital cores obtained by LSCM method is low. This leads to too little information on the obtained 3D digital core, and the results of numerical simulation using 3D digital cores are not representative. Figure from [43].


## 2.2. Imaging Methods Based on a Serial Section Rock Overlay

### 2.2.1. Ordinary Sequence Two-Dimensional (2D) Slice Superposition Method

Rock imaging methods based on a serial section rock overlay are destructive imaging techniques, of which two exist: the ordinary sequence 2D slice superposition method and the focused ion beam scanning electron microscopy (FIB-SEM) physical imaging method [46, 47]. In the former, the rock sample is prepared by polishing it to obtain a relatively flat plane, after which a high-magnification microscope is used to photograph the polished rock sample surface to obtain a microscopic image of the core surface. Subsequently, a layer of rock specimen is cut parallel to the polishing surface, and the cut rock specimen is further polished and photographed with a high-magnification microscope. This process is repeated until a 3D digital core with the desired thickness is obtained. Finally, experimental images of the core are obtained by combining these photographs. This method yields nanoscale and high-resolution core images [48–50].


### 2.2.2. FIB-SEM Physical Imaging Method

The ordinary sequence 2D slice superposition method can obtain a higher-resolution 3D digital core, but this technique uses an electron beam to polish the sample surface, which produces static electricity at the surface; therefore, imaging the surface is not beneficial. As an alternative, FIB-SEM was first developed in 1988. The FIB-SEM system utilizes two beams and is thus referred to as the dual-beam system [51]. FIB-SEM can be simply understood as a single-beam FIB combined with an SEM. Differently from the ordinary-sequence 2D slice superposition method, the use of ion-beam polishing has the advantage of producing less static electricity and better imaging quality. Therefore, the FIB-SEM imaging method is more applicable to this study than the ordinary sequence 2D slice superposition method.Although the resolution can reach the nanoscale based on the rock serial section overlay imaging method, the FIB-SEM technique is problematic because it destroys rock samples and has a slow modelling speed (Figure2). For a sample with a surface area of 50μm×50μm, the FIB abrasion of a sample with a thickness of 0.1 μm can take minutes, and only 5-20 images can be scanned per hour, which substantially slows the modelling process. In addition, the faster the sample erodes, the poorer the polishing effect becomes; thus, the quality of the image is affected. Therefore, increasing the speed of the serial section rock overlay imaging method is difficult [52, 53]. Another disadvantage of this method is the cutting and subsequent destruction of the rock sample. In addition, when a sample is recut, the thickness between adjacent slices is directly reduced during the cutting, and the pore structure of the rock sample is destroyed. Consequently, an inaccurate pore structure is obtained for the rock sample. This phenomenon results in a relatively poor pore structure after splicing that does not completely reflect the original pore structure. These three disadvantages restrict the development of this method. However, the accuracy of this technique is still better than that of focused scanning, and the nanoscale resolution is favourable for the identification of organic pores in shale rocks. Therefore, this method is used more often than focused scanning [54–57].Figure 2
The multicomponent 3D digital core reconstructed by the FIB-SEM method. Different components are reflected in different greyscales. The 3D digital core obtained by this method has a high resolution, but the disadvantage is that the scanned sample will be destroyed.


### 2.2.3. X-Ray CT Scanning Method

The basic principle of CT, which was proposed by G.N. Hounsfield in 1969 [58], involves the use of X-ray scanning to obtain the thickness of an object as the detector receives the X-rays through the sample. The X-rays are converted into electric signals that consist of visible light through a photoelectric converter and are then converted into digital signals through an analogue-digital converter; finally, the signals are processed by a computer to obtain 3D greyscale volume data. During the actual scanning process, the sample is subdivided into voxels of the same volume, and the X-ray attenuation coefficient of each voxel is determined by scanning and arrangement in a matrix, which is then converted into greyscale volume data. The attenuation coefficient is the key to distinguishing the material components of a digital core. Typically, the number of X-rays that are absorbed by a sample depends on the density of the components in the sample. Therefore, the absorption coefficient reflects the material composition of the sample.The attenuation coefficient of X-rays that pass through a component can be determined by Beer’s law as follows:(1)I=I0⋅exp−μx,where I0 is the incident X-ray intensity, I is the X-ray intensity after attenuation, μ is the linear attenuation coefficient, and x is the length of the X-ray penetration path.If the sample consists of many components (such as rocks with more complex minerals), the above equation becomes the following:(2)I=I0⋅exp∑i−μixi,where i represents the i−th component.In the early 1980s, the first microcomputed tomography system (micro-CT) was developed to increase the vertical resolution to the level of the horizontal resolution [59]. In the early 1990s, Dunsmuir proposed the application of the micro-CT method in the field of petroleum exploration [60]. Later, Coenen et al. used micro-CT techniques to construct submicron digital cores [61]. Currently, two micro-CT systems are commonly employed: a benchtop micro-CT scanning system that generates X-rays with an industrial X-ray generator and a synchronous acceleration micro-CT system that uses a synchrotron as an X-ray generator [62]. Generally, the digital core resolution that is obtained by a synchronous acceleration micro-CT scanning system is higher than that from a desktop micro-CT scanning system, but the resolution of a micro-CT scanning system is typically still within the micrometre to submicrometre range. Figure 3 shows a 3D sandstone core that was constructed from CT scanning by Arns et al. at the Australian National University, and Figure 4 shows the Fontainebleau sandstone that was constructed with CT scanning at Brookhaven National Laboratory.Figure 3
The use of X-ray CT scanning method is the establishment of a 3D digital core. The transparent part is the skeleton; the other part is the pores.Figure 4
The Fontainebleau sandstone dataset with a binary segmentation. The red part is the skeleton; the black part is the pores. The total size of the dataset is288×288×300 and a voxel edge length is 7.5 μm. Figure from [63].In theory, the X-ray resolution limit should be 0.005 nm, but the original micro-CT technique cannot probe at a resolution greater than the submicron scale. The key restriction is the size of the X-ray source and the pixel size of the detector. The detector size of the X-ray source cannot reach the nanometre scale. Fortunately, the resolution of a current radiation source can reach the nanometre scale, and the size of the radiation source can be constrained to achieve the resolution of nanometre-scale CT (nano-CT) [64, 65]. After the appearance of the nano-CT technique, nanometre CT systems based on visible light optical systems, synchrotron radiation sources, and X-ray optical systems have been proposed. For example, Lawrence Berkeley National Laboratory obtained 15 nm resolution CT images. With continuous improvements in CT scanning technology, the resolution of this method can be further improved.Compared to the above two methods, the greatest advantage of the X-ray CT scanning method is that the sample is not damaged during the construction of the digital core. In addition, current nano-CT systems can achieve nanoscale resolution with the same accuracy as FIB-SEM systems. Table1 shows the sample size, instrument resolution, and main application areas of the CT and FIB-SEM analysis techniques.Table 1
Instrument parameters and applications for 3D digital core reconstruction.
Instrument type
Sample size
Instrument resolution
Main application areas
Micro CT
Core column scanning(diameter 15-25 mm)
12-30μm
Geological description, heterogeneity analysis
Core column scanning(diameter 1-25 mm)
0.5-12μm
Fine description, pore network generation
Nano-CT
Less than 1 mm
65-150 nm
For the scanning of complex structural samples such as shale, microcrystal, tight gas rock, and analysis of pore network
FIB-SEM
Less than 1 mm
2-150 nmAmong the abovementioned methods, the X-ray CT scanning technique is currently the most popular. However, this approach also suffers from some problems. For instance, the sample size of a CT experiment is too small. For rocks with a high degree of heterogeneity, the results from CT scanning are not representative, so one must reconstruct the multiscale digital core. Nevertheless, the applicability of CT experiments is generally broader than that of FIB-SEM experiments.


### 2.2.4. Magnetic Resonance Imaging (MRI) Method

The magnetic resonance imaging (MRI) method was first proposed in 1973 by Lauterbur [66]. Subsequently, MRI methods were gradually applied in the medical and biological communities, among others. In 1995, Xiao applied the MRI method to the field of petroleum engineering for the first time [67], after which this technique was widely employed in both petrophysics and petroleum engineering [68–73]. The use of the MRI method to image fluids in rocks to reflect their structures has constituted an important guide in the visual determination of microscopic pore structures.A nucleus with a nonzero magnetic moment undergoes a phenomenon known as NMR; that is, Zeeman splitting occurs in the spin-energy level of the external magnetic field, and the resonance absorbs radiation that is associated with radio waves at a certain frequency. NMR signals are generated during NMR within the nucleus. These NMR signals are collected by a special device, and a 2D inverse Fourier transform is used to obtain a reconstructed image of a spatial point. The entire target object can be imaged by collecting the signals at each point in the range of the target object.Compared to the LSCM, FIB-SEM, and X-ray CT methods, the MRI method is limited by a very obvious disadvantage: an excessively low imaging resolution [74, 75]. For tight sandstones and shale rocks, obtaining a large amount of microscopic information with the MRI method is difficult. In addition, NMR experiments require that the pores of the rock be saturated with fluids, which is inherently difficult to achieve for very dense rocks. Moreover, the type of fluid within the rock can affect the imaging results [76, 77]. The NMR inversion algorithm is another factor that influences the imaging quality; accordingly, research on NMR inversion algorithms is ongoing [78–80]. In summary, MRI is not currently employed to construct digital cores because of its limitations.Many relatively simple pore network models are also available [81, 82], and experiments such as mercury intrusion and one-dimensional NMR can be used to image the porosity of a rock [83, 84]. However, the former cannot accurately reflect the heterogeneity of rocks. Moreover, the pore shape and overall pore network in a simple pore network model are fixed, and the pores are arranged in a given orientation; these principles are in contrast to the idea that the digital core must reflect the rock features, including the heterogeneity, to the greatest extent possible. Unfortunately, experiments such as mercury intrusion reflect the pore structure only from the side of the sample, which is not intuitive and has multiple solutions. Therefore, the above two methods cannot be called digital core construction methods.Overall, the most promising and suitable physical experimental method that is currently available for the reconstruction of a digital core is CT scanning. The FIB-SEM reconstruction method can also obtain digital cores with relatively high resolution, but the samples will be destroyed in the process. Nevertheless, as the resolution of these instruments continues to increase, the accuracy of physical experimental techniques in the reconstruction of 3D digital cores will also continue to improve.


## 2.2.1. Ordinary Sequence Two-Dimensional (2D) Slice Superposition Method

Rock imaging methods based on a serial section rock overlay are destructive imaging techniques, of which two exist: the ordinary sequence 2D slice superposition method and the focused ion beam scanning electron microscopy (FIB-SEM) physical imaging method [46, 47]. In the former, the rock sample is prepared by polishing it to obtain a relatively flat plane, after which a high-magnification microscope is used to photograph the polished rock sample surface to obtain a microscopic image of the core surface. Subsequently, a layer of rock specimen is cut parallel to the polishing surface, and the cut rock specimen is further polished and photographed with a high-magnification microscope. This process is repeated until a 3D digital core with the desired thickness is obtained. Finally, experimental images of the core are obtained by combining these photographs. This method yields nanoscale and high-resolution core images [48–50].


## 2.2.2. FIB-SEM Physical Imaging Method

The ordinary sequence 2D slice superposition method can obtain a higher-resolution 3D digital core, but this technique uses an electron beam to polish the sample surface, which produces static electricity at the surface; therefore, imaging the surface is not beneficial. As an alternative, FIB-SEM was first developed in 1988. The FIB-SEM system utilizes two beams and is thus referred to as the dual-beam system [51]. FIB-SEM can be simply understood as a single-beam FIB combined with an SEM. Differently from the ordinary-sequence 2D slice superposition method, the use of ion-beam polishing has the advantage of producing less static electricity and better imaging quality. Therefore, the FIB-SEM imaging method is more applicable to this study than the ordinary sequence 2D slice superposition method.Although the resolution can reach the nanoscale based on the rock serial section overlay imaging method, the FIB-SEM technique is problematic because it destroys rock samples and has a slow modelling speed (Figure2). For a sample with a surface area of 50μm×50μm, the FIB abrasion of a sample with a thickness of 0.1 μm can take minutes, and only 5-20 images can be scanned per hour, which substantially slows the modelling process. In addition, the faster the sample erodes, the poorer the polishing effect becomes; thus, the quality of the image is affected. Therefore, increasing the speed of the serial section rock overlay imaging method is difficult [52, 53]. Another disadvantage of this method is the cutting and subsequent destruction of the rock sample. In addition, when a sample is recut, the thickness between adjacent slices is directly reduced during the cutting, and the pore structure of the rock sample is destroyed. Consequently, an inaccurate pore structure is obtained for the rock sample. This phenomenon results in a relatively poor pore structure after splicing that does not completely reflect the original pore structure. These three disadvantages restrict the development of this method. However, the accuracy of this technique is still better than that of focused scanning, and the nanoscale resolution is favourable for the identification of organic pores in shale rocks. Therefore, this method is used more often than focused scanning [54–57].Figure 2
The multicomponent 3D digital core reconstructed by the FIB-SEM method. Different components are reflected in different greyscales. The 3D digital core obtained by this method has a high resolution, but the disadvantage is that the scanned sample will be destroyed.


## 2.2.3. X-Ray CT Scanning Method

The basic principle of CT, which was proposed by G.N. Hounsfield in 1969 [58], involves the use of X-ray scanning to obtain the thickness of an object as the detector receives the X-rays through the sample. The X-rays are converted into electric signals that consist of visible light through a photoelectric converter and are then converted into digital signals through an analogue-digital converter; finally, the signals are processed by a computer to obtain 3D greyscale volume data. During the actual scanning process, the sample is subdivided into voxels of the same volume, and the X-ray attenuation coefficient of each voxel is determined by scanning and arrangement in a matrix, which is then converted into greyscale volume data. The attenuation coefficient is the key to distinguishing the material components of a digital core. Typically, the number of X-rays that are absorbed by a sample depends on the density of the components in the sample. Therefore, the absorption coefficient reflects the material composition of the sample.The attenuation coefficient of X-rays that pass through a component can be determined by Beer’s law as follows:(1)I=I0⋅exp−μx,where I0 is the incident X-ray intensity, I is the X-ray intensity after attenuation, μ is the linear attenuation coefficient, and x is the length of the X-ray penetration path.If the sample consists of many components (such as rocks with more complex minerals), the above equation becomes the following:(2)I=I0⋅exp∑i−μixi,where i represents the i−th component.In the early 1980s, the first microcomputed tomography system (micro-CT) was developed to increase the vertical resolution to the level of the horizontal resolution [59]. In the early 1990s, Dunsmuir proposed the application of the micro-CT method in the field of petroleum exploration [60]. Later, Coenen et al. used micro-CT techniques to construct submicron digital cores [61]. Currently, two micro-CT systems are commonly employed: a benchtop micro-CT scanning system that generates X-rays with an industrial X-ray generator and a synchronous acceleration micro-CT system that uses a synchrotron as an X-ray generator [62]. Generally, the digital core resolution that is obtained by a synchronous acceleration micro-CT scanning system is higher than that from a desktop micro-CT scanning system, but the resolution of a micro-CT scanning system is typically still within the micrometre to submicrometre range. Figure 3 shows a 3D sandstone core that was constructed from CT scanning by Arns et al. at the Australian National University, and Figure 4 shows the Fontainebleau sandstone that was constructed with CT scanning at Brookhaven National Laboratory.Figure 3
The use of X-ray CT scanning method is the establishment of a 3D digital core. The transparent part is the skeleton; the other part is the pores.Figure 4
The Fontainebleau sandstone dataset with a binary segmentation. The red part is the skeleton; the black part is the pores. The total size of the dataset is288×288×300 and a voxel edge length is 7.5 μm. Figure from [63].In theory, the X-ray resolution limit should be 0.005 nm, but the original micro-CT technique cannot probe at a resolution greater than the submicron scale. The key restriction is the size of the X-ray source and the pixel size of the detector. The detector size of the X-ray source cannot reach the nanometre scale. Fortunately, the resolution of a current radiation source can reach the nanometre scale, and the size of the radiation source can be constrained to achieve the resolution of nanometre-scale CT (nano-CT) [64, 65]. After the appearance of the nano-CT technique, nanometre CT systems based on visible light optical systems, synchrotron radiation sources, and X-ray optical systems have been proposed. For example, Lawrence Berkeley National Laboratory obtained 15 nm resolution CT images. With continuous improvements in CT scanning technology, the resolution of this method can be further improved.Compared to the above two methods, the greatest advantage of the X-ray CT scanning method is that the sample is not damaged during the construction of the digital core. In addition, current nano-CT systems can achieve nanoscale resolution with the same accuracy as FIB-SEM systems. Table1 shows the sample size, instrument resolution, and main application areas of the CT and FIB-SEM analysis techniques.Table 1
Instrument parameters and applications for 3D digital core reconstruction.
Instrument type
Sample size
Instrument resolution
Main application areas
Micro CT
Core column scanning(diameter 15-25 mm)
12-30μm
Geological description, heterogeneity analysis
Core column scanning(diameter 1-25 mm)
0.5-12μm
Fine description, pore network generation
Nano-CT
Less than 1 mm
65-150 nm
For the scanning of complex structural samples such as shale, microcrystal, tight gas rock, and analysis of pore network
FIB-SEM
Less than 1 mm
2-150 nmAmong the abovementioned methods, the X-ray CT scanning technique is currently the most popular. However, this approach also suffers from some problems. For instance, the sample size of a CT experiment is too small. For rocks with a high degree of heterogeneity, the results from CT scanning are not representative, so one must reconstruct the multiscale digital core. Nevertheless, the applicability of CT experiments is generally broader than that of FIB-SEM experiments.


## 2.2.4. Magnetic Resonance Imaging (MRI) Method

The magnetic resonance imaging (MRI) method was first proposed in 1973 by Lauterbur [66]. Subsequently, MRI methods were gradually applied in the medical and biological communities, among others. In 1995, Xiao applied the MRI method to the field of petroleum engineering for the first time [67], after which this technique was widely employed in both petrophysics and petroleum engineering [68–73]. The use of the MRI method to image fluids in rocks to reflect their structures has constituted an important guide in the visual determination of microscopic pore structures.A nucleus with a nonzero magnetic moment undergoes a phenomenon known as NMR; that is, Zeeman splitting occurs in the spin-energy level of the external magnetic field, and the resonance absorbs radiation that is associated with radio waves at a certain frequency. NMR signals are generated during NMR within the nucleus. These NMR signals are collected by a special device, and a 2D inverse Fourier transform is used to obtain a reconstructed image of a spatial point. The entire target object can be imaged by collecting the signals at each point in the range of the target object.Compared to the LSCM, FIB-SEM, and X-ray CT methods, the MRI method is limited by a very obvious disadvantage: an excessively low imaging resolution [74, 75]. For tight sandstones and shale rocks, obtaining a large amount of microscopic information with the MRI method is difficult. In addition, NMR experiments require that the pores of the rock be saturated with fluids, which is inherently difficult to achieve for very dense rocks. Moreover, the type of fluid within the rock can affect the imaging results [76, 77]. The NMR inversion algorithm is another factor that influences the imaging quality; accordingly, research on NMR inversion algorithms is ongoing [78–80]. In summary, MRI is not currently employed to construct digital cores because of its limitations.Many relatively simple pore network models are also available [81, 82], and experiments such as mercury intrusion and one-dimensional NMR can be used to image the porosity of a rock [83, 84]. However, the former cannot accurately reflect the heterogeneity of rocks. Moreover, the pore shape and overall pore network in a simple pore network model are fixed, and the pores are arranged in a given orientation; these principles are in contrast to the idea that the digital core must reflect the rock features, including the heterogeneity, to the greatest extent possible. Unfortunately, experiments such as mercury intrusion reflect the pore structure only from the side of the sample, which is not intuitive and has multiple solutions. Therefore, the above two methods cannot be called digital core construction methods.Overall, the most promising and suitable physical experimental method that is currently available for the reconstruction of a digital core is CT scanning. The FIB-SEM reconstruction method can also obtain digital cores with relatively high resolution, but the samples will be destroyed in the process. Nevertheless, as the resolution of these instruments continues to increase, the accuracy of physical experimental techniques in the reconstruction of 3D digital cores will also continue to improve.


## 3. Status of 3D Digital Core Technology Stochastic Reconstruction by 2D Slices

Digital cores can be acquired with sufficient accuracy by 3D digital core reconstruction methods based on a physical experimental method. However, constructing large-scale digital cores with physical experimental reconstruction methods, which greatly influence the accuracy of simulation experiments on rocks with a high degree of heterogeneity, remains difficult. More importantly, physical experiments are expensive and time-consuming [85]. Unfortunately, the conventional abstract capillary bundle and random pore network models offer only abstract expressions of rocks and cannot be utilized to characterize their randomness. Therefore, the collection of single or multiple 2D slices through experimental methods such as SEM and thin slicing is valuable for the random reconstruction of 3D digital cores. Accordingly, the benefits of these methods have motivated many scholars to conduct in-depth and meticulous research on their applications. Existing reconstruction methods constrain and reconstruct 3D digital cores by using 2D slice morphology information and 2D slice statistics of particle size characteristics and simulating rock generation processes, among others.


### 3.1. Stochastic Reconstruction Based on Statistics

#### 3.1.1. Stochastic Reconstruction Based on 2D Slice Statistical Constraint Features

The stochastic reconstruction method is the most common method for reconstructing 3D digital cores based on one or more slices of 2D slice data, various mathematical-feature statistics, and functional information, such as porosity data. Similar to facial reconstruction with abstract facial features (e.g., eye size), the stochastic reconstruction method uses the abstract mathematical features of extracted 2D slices to build the most accurate 3D digital cores [86–90].(1) Gaussian Field Method. In 1974, Joshi was the first to use stochastic information that was derived from the statistics of 2D slices to randomly reconstruct rocks [91]. Joshi proposed a reconstruction technique known as the Gaussian field method, which can be roughly divided into three steps. First, a Gaussian field is randomly generated, and the Gaussian field is filled with independent Gaussian variables. Then, the variables are linearly transformed to correlate the variables. The constraints that are used during the transformation are the core porosity values and the two-point correlation function. Finally, the Gaussian field from the second step is transformed into a digital core by using a nonlinear transformation. However, Joshi built only 2D digital core models at the time.In 1984, Quiblier [92] improved upon the primitive Gaussian field-reconstruction method and constructed 3D digital cores. Subsequently, Adler et al. [93] introduced periodic boundary conditions and the Fourier transform to improve the speed and efficiency of the model, respectively. Furthermore, Hilfer [94] added porosity and probability distribution functions to the model to improve its accuracy. Torquato and Lu [95] used the chord-length distribution function to describe pore space features. Ioannidis et al. [96] added the fast Fourier transform to improve the microscopic pore structure accuracy of the constructed digital core. Giona and Adrover [97] used linear filters to improve the accuracy of the reconstructed model. Roberts [98] replaced the two-point probability function with a multipoint probability function. Most recently, Liang et al. [99] improved the Fourier transform by using a truncated Gaussian method.The abovementioned methods clearly improved the use of the Gaussian field reconstruction method; however, this technique has been criticized because of the poor connectivity within the constructed digital core. The connectivity of the reconstructed 3D digital core does not greatly improve with the use of any of the abovementioned methods; thus, further development of the Gaussian field reconstruction method is necessary. Generally, the advantage of the Gaussian field method is its fast modelling speed. If the accuracy of the digital core is not sufficient, then no quantitative calculation is required if the core to be reconstructed is a high-porosity rock. If a more accurate digital core is required, then the number of constraints must be increased; that is, the evaluation function must be enhanced. However, if additional constraints are included, the modelling speed of the Gaussian field method diminishes, so long-distance connectivity problems cannot be solved, and the accuracy becomes insufficient. To obtain a more accurate 3D digital core in this scenario, one should use a digital 3D core reconstruction method based on an optimization algorithm (such as the simulated annealing (SA) algorithm) with similar principles. Accordingly, the Gaussian field reconstruction approach may be replaced eventually.(2) Stochastic Search Algorithm Based on a Majority Operator. In 2007, Zhao et al. proposed a random search method based on a majority operator and evaluated the permeability of the reconstructed digital core [100]. The main modelling process of the majority operator is divided into four steps. First, the 2D slice two-point probability and linear path functions are calculated. Second, datasets that are the same size as the 3D digital cores are randomly generated (for example, a 300×300×300 voxel digital core is randomly generated when reconstructing a 300×300×300 voxel digital core), and the difference between the two-point probability and linear path functions is calculated as the objective function value. Third, the randomly generated 3D digital cores are optimized based on the majority operator-based random search strategy. Seven main models are used for the optimizing operators: D2Q9, D3Q13, D3Q15, D3Q19, D1Q3, D2Q5, and D3Q27 [101]. Assuming that the D3Q27 model is selected, a voxel that points to a hole can be randomly selected in a digital core; then, the binary states of its 26 neighbourhoods (i.e., the porosity or rock skeleton) can be judged. If a point represents a pore among the 26 neighbourhoods that exceeds a certain value (usually set as a constant), the point is selected again; otherwise, the point is selected and replaced with a skeleton point. After each replacement is completed, the target function value is recalculated until the target function value falls below a given value or the target function value no longer changes over several iterations.Only scarce literature is available regarding the multiple-operator random reconstruction method, which is more suitable for the construction of 3D digital cores with large pores (usually more than 20% porosity) [102] because this approach suffers from many weaknesses. First, this method is optimized with a fixed law, which is very slow, and thus easily falls into extreme values, making it difficult for the model to continue the optimization (i.e., SA and other metaheuristics are random and are therefore the best approach to solve the traditional optimization problem). Second, the adjustment of this method is partial, not integral. This partial adjustment can guarantee only the details, while the overall long-range connectivity is difficult to guarantee. Although the majority operator method works better than the Gaussian field method, it still does not address the holistic, random, and long-range connectivity that is required for optimized digital core-reconstruction methods. However, the multiple-operator random reconstruction method can be used as an optimization solution for reconstructing digital cores (described in Section 3.3).(3) Simulated Annealing (SA) Method. In 1997, Hazlett [103] first proposed a new digital core-reconstruction method based on the SA optimization algorithm. After evaluating the 2D slice porosity, the two-point probability function, and the linear path, a 3D digital core that has the same porosity as the 2D slice of the rock is first constructed, after which the digital core is optimized by using an SA optimization algorithm. During each iteration of the optimization, a pore voxel point and skeleton voxel point are randomly selected and exchanged, and the objective function is calculated. If the value of the function decreases, the 3D digital core is updated, the iteration is continued until it terminates, and a reconstructed 3D digital core is obtained. Figure 5 shows the flow chart for reconstructing a 3D digital core with the SA method. Figure 6 shows a 3D digital sandstone core that was obtained with the SA optimization algorithm in Figure 6, in which blue colors denote pores and red colors represent the skeleton.Figure 5
SA method reconstruction of 3D digital core reconstruction.Figure 6
Reconstruction of 3D digital cores by the SA method. The blue part is the pore; the red part is the skeleton. It can be seen that the constructed 3D digital core has better pore connectivity and is similar to the real pore space. Figure reprinted from [104].Subsequently, Yeong, Torquato and Manwart [105–107] proved that the SA method can approximate the evaluation function. Eschricht et al. [108] improved the SA method and the accuracy of the reconstructed model. After building digital cores with the SA method, Zhao and Yao [109] suggested that good evaluation parameters can improve the reconstruction accuracy and consequently proposed an improved method to speed up the model optimization process. Considering the disadvantages of the SA optimization method, Teng et al. [110] proposed a 3D image reconstruction method based on particle swarm optimization (PSO); an improved reconstruction was obtained compared to the Gaussian field reconstruction method and the SA method. Ju et al. [111] proposed an improved SA method that enhanced the modelling effect. Likewise, Zou et al. [112] improved the SA method and proposed a very fast SA optimization method, which was applied to reconstruct shale reservoirs. Frączek et al. [113] proposed a multiscale SA reconstruction method, and Zeng et al. used the SA reconstruction method to reconstruct gas-hydrate digital cores, achieving very good results. Lin et al. [114] proposed a multiscale SA reconstruction method and reconstructed digital carbonate-rock cores with good results. In 2018, Capek proposed an improved SA method to further improve the reconstruction accuracy [115].The effectiveness of reconstructions from the SA method mainly depends on the efficiency of the algorithm and the use of an appropriate evaluation function. With the same evaluation function, the effect of the construction of the SA method is better than that of the Gaussian field method as well as that of the majority operator. The key issue at present is that strong global optimization is the most important component to improve the long-range connectivity regardless of the optimization model. Thus, a large number of evaluation functions that reflect the different details of digital cores can be added to improve the prediction effect. Among the several methods that are currently employed for this task, metaheuristic optimization algorithms, including SA optimization algorithms, are the most promising because of their high optimization efficiency and certain randomness (thereby avoiding extreme values and improving the final optimization effect). In Section4, we also mention that the prediction may be further improved if a new metaheuristic optimization algorithm is used. In addition, the combination of multiple optimization algorithms based on their unique characteristics can further improve the digital core-reconstruction effect.(4) Discussion of the Statistically Constrained Stochastic Reconstruction of Mathematical Features. Many methods are based on the statistically constrained stochastic reconstruction of mathematical features. In this paper, four common methods are introduced. Among them, the most studied is the SA method, which is a reconstruction method that we consider to be very promising. The digital cores that are reconstructed by the SA method boast better connectivity than reconstructed 3D digital cores from the stochastic search algorithm based on the majority operator and Gaussian field methods; thus, the constraint function of the evaluation function can be better executed. In this section, we focus on a discussion of SA methods. To reconstruct a digital core, the effect of an algorithm depends on the required time for the reconstruction and the approximation degree (error size) of the evaluation function. To improve these two components, an appropriate initial model and the algorithm itself should be reasonably improved. The appropriate initial model can be based on a combination of different methods, which are discussed in detail in Section 3.3.The SA algorithm, which was first proposed in 1953, is a classic optimization algorithm that regards the accuracy of the 3D digital core reconstruction as an optimization problem. SA algorithm optimization works very well for rocks with high porosity and high permeability and rocks with uncomplicated pore structures. However, the speed and accuracy are greatly reduced when reconstructing large-scale rocks with complex pore structures [116], indicating that the classical SA algorithm cannot meet complex optimization requirements. Many scholars that have engaged in 3D digital core-reconstruction research have not addressed any improvements to this algorithm itself. Instead, scholars who have focused on optimization research have made improvements to various individual aspects of the SA algorithm, such as the ease with which the algorithm falls into local minima [117, 118], its long execution time [119–125], and the sensitivity of its parameters [126–128]. Since the concept of the metaheuristic optimization algorithm [103] was proposed, research progress regarding this algorithm has been very rapid. Numerous new metaheuristic optimization algorithms based on various phenomena have been proposed. Compared to these metaheuristics, the SA algorithm is primitive. Table 2 summarizes some of the classic metaheuristic optimization algorithms and the latest metaheuristic optimization algorithms. Figure 7 shows a 3D digital core that was reconstructed by PSO and the SA algorithm with the same evaluation function.Table 2
Some of the metaheuristic optimization algorithms.
Imitated object
Representative algorithms
The evolution of nature
Genetic algorithm [129], differential evolution algorithm [130]
Human behaviour
Immune algorithm [131], brainstorm algorithm [132], fireworks algorithm [133]
Physical characteristics
Simulated annealing algorithm [134], intelligent water drop algorithm [135], binary black hole algorithm [136]
Animal behaviour
Particle swarm optimization algorithm [137], ant algorithm [138], bat algorithm [139], krill herd algorithm [140], whale optimization algorithm [141]
Plant behaviour
Invasive weed algorithm [142], rain forest algorithm [143], root growth model [144]Figure 7
Particle swarm optimization (a), Fourier transform-based Gaussian field method (b), and simulated annealing method (c) reconstruction results of 3D digital core. Among them, the 3D digital core with the same background color is the skeleton, and the white part is the pores. It can be observed that the reconstruction effect of particle swarm optimization is better than that of the simulated annealing algorithm, and the obtained 3D digital core has better connectivity in the pore space. Figure from [109].
(a)
(b)
(c)Table2 demonstrates the evident diversity of heuristic optimization algorithms; moreover, only a small number of the many available optimization algorithms are shown in Table 2. Hundreds of optimization algorithms can be used, even if the improved algorithm is not calculated; in addition, no improved algorithms based on the abovementioned metaheuristic optimization algorithm are included in this table. Many of these optimization algorithms have been proven superior to the 35-year-old SA method in terms of the search speed and optimization effect. For the reconstruction of 3D digital cores, one should use a variety of the most recently developed optimization algorithms to determine which is the most suitable for digital core-reconstruction problems and research and to improve the search ability of each algorithm as well as its optimization speed in order to satisfy the needs of higher evaluation functions. This approach is also the best to improve the accuracy of the stochastic reconstruction method based on the statistical constraint of mathematical features. An interdisciplinary method should be used to improve the reliability of a reconstructed 3D digital core. Alternatively, we can combine the characteristics of multiple optimization methods to jointly reconstruct a digital core. For example, the Gaussian field method is used first to optimize the initial stochastic model; at this time, fewer evaluation functions are selected. The results are then transferred to the optimization algorithm, and more evaluation functions are added to refine the digital core modelling process, thereby increasing the long-range connectivity and reliability of the digital core. Finally, a multituning method is used to further fine-tune the model, and then the model more closely resembles the actual rock. This approach can also be regarded as a hybrid optimization method. In summary, the advantage of the 2D slice method based on the statistical constraint of mathematical features is that it considers the integrity of the digital core during its reconstruction; in addition, this technique first generates the 3D random data and then performs the optimization. However, the disadvantages of this approach are obvious. Morphological characteristics such as sex are obviously insufficient and depend on mathematical algorithms.


### 3.2. Stochastic Reconstruction Method Based on the Statistical Constraints of 2D Slice Morphological Characteristics

A stochastic reconstruction method based on the statistical constraints of 2D slice morphological characteristics can be used to calculate the morphological features of 2D slices (specifically, the relationships between the positions and relative positions of voxels). This method can also be used for 3D reconstruction purposes based on the morphological features of 2D slices, but it differs from the stochastic method based on the statistical constraints of mathematical features, which considers the correlations between voxels with regard to the connectivity and morphological features of pores. Representative reconstruction algorithms include the Markov chain Monte Carlo (MCMC) reconstruction method and the stochastic reconstruction method based on the statistical constraints of slice morphology characteristics; these algorithms are used to calculate the morphological features of 2D slices (specifically, the relationships between the positions and relative positions of voxels). A reconstructed 3D digital core is then built based on the morphological features of these 2D slices. This approach is different from the stochastic method based on statistically constrained mathematical features, which considers the correlations between voxels with regard to the connectivity and morphological features of pores. Representative reconstruction algorithms include the MCMC reconstruction method and multipoint geostatistics (MPS) [145].


#### 3.2.1. Sequence Indicator Simulation Method

In 2003, Keehm [146] proposed a new reconstruction method, namely, a 3D digital core stochastic reconstruction method based on the sequence indicator simulation method from geostatistics. Sequence indicator simulation is a method that indicates the probability distribution field of data by using a directed kriging interpolation technique in combination with a conditional stochastic simulation. In this method, the variogram and porosity values of the 2D slice are first calculated. Then, the path of each network node is randomly assessed, and kriging interpolation is applied to the indicator variable to determine the probability that this variable is a discrete variable at the node. Then, the discrete variable order is determined, and a random number is generated, which determines the random variable type at that point. The nodes are repeatedly updated until each node is simulated, and a 3D digital core is ultimately obtained. In 2007, Zhu and Tao [147] reconstructed a 3D digital core by using cast sheet data. The variation function of the reconstructed digital core is considered to correspond well with the variation function that is calculated from the reconstructed casting slice. The sequence indicator simulation method produces reconstructed 3D digital cores that are very similar to real digital cores, so this method can be applied to the simulation of electrical and elastic rock properties. In 2011, Liu and Mu [148] proposed an improved sequence indicator simulation method and applied it to calculate the permeability of a digital core.Unfortunately, the sequence indicator simulation method does not fundamentally solve the reconstructed digital core connectivity problem. The greatest advantage of the sequence indicator simulation method is that it has certain integrity and correlations between voxels, which is very important for the reconstruction of digital cores to ensure long-range connectivity. The relevance of each voxel is greatly considered. In contrast, the weakness of the sequence indicator simulation method is that the rules to establish the digital core are directly fixed from the image, which causes two hidden problems. First, when using this method for the reconstruction of a digital core, the reconstruction accuracy is directly derived from the image. For images with very high heterogeneity, the image representation is weak, and the reconstruction effect is influenced accordingly. Second, because the evaluation function is relatively fixed, developing it further is difficult, so the model has weaker potential. The above problem is a common problem of modelling methods that are based on the morphological features of images, but the disadvantages of the sequence indicator simulation method are more serious. Compared to the MCMC and MPS methods, the digital cores that are reconstructed by the sequence indicator simulation method are also less reliable and are especially weaker overall than those that are reconstructed by the MPS method, which considers several consecutive voxel points in the vicinity and thereby more fully considers the digital core. The subsequent implementation of the simulation method is significantly worse than the implementation of both the MCMC and MPS methods. In summary, the sequence indicator simulation method does not fundamentally solve the reconstructed digital core connectivity problem.


#### 3.2.2. Markov Chain Monte Carlo (MCMC) Method

In 2004, the MCMC reconstruction method was proposed by Wu et al., who applied this technique to the 2D reconstruction of soil structure [149]. In 2006, Wu et al. extended the 2D digital core-reconstruction method into three dimensions and proposed an MCMC-based 3D digital core-reconstruction method [150]. This technique introduces the Markov chain and traverses all two-point and five-point neighbourhoods in a 2D slice image to calculate the conditional probability of the neighbourhood template; this method also determines the state of each point in the reconstructed image based on conditional probabilities. Markov chains describe a sequence of states (analogous to the notion of a template). The state value of each position in the sequence depends on the state of the previous position, and the probability of describing this state is called the transition probability [151, 152]. This reconstruction method can be roughly divided into four steps. First, we establish a Markov chain with a stable probability distribution function by means of traversal scanning. When the chain is performed over a sufficient distance in the image, an important statistical feature of the original image is obtained. After the conditional probability is obtained, the image reconstruction can be performed. Then, the first row of the first layer is reconstructed, and the average porosity is determined by using 2D slices in the x, y, and z directions (slices can be reconstructed in three directions from a single 2D slice). From the second voxel, a two-point neighbourhood template is used to reconstruct the three-point neighbourhood system from the third point to reconstruct the first row of voxel values. The third step is to reconstruct the first layer of the 3D digital core. The first point is reconstructed with two-point neighbourhood templates, and the second point is reconstructed by four-point neighbourhood templates. The remaining points are reconstructed by using the five- and six-point neighbourhood systems. Subsequently, the remaining layers of the 3D digital core are reconstructed, and the first layer is reconstructed vertically by using the previous steps. Then, the first voxel of the second layer utilizes seven- and eight-point neighbourhood systems for reconstruction. Finally, the second line of the remaining components of the 3D digital core can be reconstructed with 10- and 11-point neighbourhood systems. A diagram of the MCMC reconstruction method is shown in Figure 8, and Figure 9 displays a digital core of the Berea sandstone that was reconstructed with the MCMC method.Figure 8
Reconstruction process of 3D digital cores using MCMC reconstruction method. (a) Shows reconstruction the first voxel, (b) shows reconstruction the first row, (c) shows reconstruction the first level, and (d) shows reconstruction the remaining layers of the 3D digital core. Figure from [153].
(a)
(b)
(c)
(d)Figure 9
Reconstruction of 3D digital cores using MCMC reconstruction method. Among them, black is the pore part, and white is the skeleton part. The 3D digital cores obtained by the MCMC method have good connectivity of pore space and more concentrated pore space.Wang et al. [154] used the MCMC method to reconstruct the pores of shale rocks in two dimensions. Nie et al. [154] modelled a shale-gas reservoir with a large-scale 3D digital core by using the 3D MCMC method. Guo et al. [155] used a random walk algorithm to simulate the NMR response of a sandstone reservoir in a reconstructed 3D digital core based on the MCMC algorithm. Overall, the MCMC-reconstructed sandstone digital cores match the actual cores, and the NMR response values are similar.The MCMC method is a very reliable digital core-reconstruction technique. The constructed digital core boasts good connectivity. However, MCMC-reconstructed digital cores have weak anisotropy, and the pore throat radius distribution is very concentrated. These factors are the major disadvantages for reservoirs with strong heterogeneity. Compared to the sequential indicator simulation method, the MCMC method considers the voxel configuration around each voxel to be determined and thus obtains a digital core with better connectivity, especially with regard to the distribution of pores, when the core is more uniform. However, the most fundamental cause of the weak anisotropy in reconstructed digital cores is that this hypothesis does not conform to a heterogeneous rock. The assumption of this technique is that the state of any point depends only on the states of a few neighbouring points. Then, after counting the probability function of the training image, the voxels that actually must be determined are determined only by points that are close to each other, whereas the distances to voxel points that are far from the points to be determined are ignored, resulting in the absence of long-distance differences and eventually producing weak anisotropy. Nevertheless, the MCMC method is generally more accurate than the sequential indicator simulation method; consequently, more studies have been performed based on the reconstruction of digital cores with the MCMC method. In conclusion, 3D digital cores that are reconstructed by the MCMC algorithm match real digital cores, but whether these cores are applicable to rocks with very complicated pore structures remains uncertain.


#### 3.2.3. Multipoint Geostatistics (MPS) Method

In 2004, Okabe and Blunt proposed a digital core-reconstruction method based on the statistical method of MPS [42]. The MPS method extracts multipoint statistical information from training images (i.e., 2D slices) and reconstructs the template and reconstruction mode; then, this mode is used to perform a 3D digital core reconstruction. Compared to the MCMC method, this method is more biased toward the reconstruction of 3D digital cores with the morphological features of 2D slices. The reconstruction of a digital core with the MPS method is divided into four steps. First, a 2D slice search template and a search tree are established. Second, the raw data are reloaded to the nearest simulation grid node and fixed during the simulation. Next, we define a path that randomly accesses all voxels. We use the search template to define conditional data events, compute the conditional probability distribution function (CPDF) based on the search tree, calculate the simulated values with the function, and generate new 2D images through iteration. Finally, we employ the reconstructed image as a new training image to generate the next layer of the image, which is then used to generate a new 3D digital core. Figure 10 shows a reconfigured flowchart of the MPS reconstruction method, and Figure 11 shows a 3D digital shale-rock core that was reconstructed with the MPS method.Figure 10
MPS reconstruction method flow chart.Figure 11
Reconstruction result of the 3D digital core using MPS method. Among them, the blue part is the pore, and the transparent part is the skeleton. It can be seen that even for shale rock, the long-range connectivity of the pore space of the 3D digital core obtained by the MPS method is still good. This shows that the MPS method can reflect the training image well and obtain the 3D digital core with long-range connectivity. Figure reprinted from [155].In 2004, Daı̈an et al. [156] used MPS to reconstruct digital cores and calculate their porosity and permeability, which matched the actual experimental results. Because the MPS method suffers from a slow reconstruction speed and occupies more memory, Wu et al. [157] proposed a solution to replace this pixel-wise distance calculation with a filter score comparison, improving the reconstruction speed by a factor of 10. In 2009, Wang et al. [158] reconstructed 3D digital cores by using 2D CT images and tested the permeability of the reconstructed 3D digital cores by using the lattice Boltzmann method. Teng et al. [159] proposed that higher multipoint statistical stationarity and lower-scale stationarity are likely to produce results that more closely resemble the real 3D structure. In 2012, Tahmasebi et al. [160] proposed a multipoint statistical method based on cross-correlation function geostatistical methods (CCSIM) and conducted continuous research to improve the method [161–165]. Straubhaar et al. [166] proposed a parallel multipoint geostatistical method based on an improved search tree structure, effectively reducing the memory footprint and increasing the computational speed of the method. In 2013, Du and Zhang [167] used multigrid templates to improve the reconstruction effect and achieved good results. In 2014, Zhang et al. [168] proposed an MPS method based on GPU acceleration. In the same year, Gao et al. [169] proposed an MPS method that uses three-step sampling. Zhang et al. [170] proposed an Isomap-based MPS reconstruction method that uses the Isomap method to reduce the dimensions and redundancy. Liu et al. [171] proposed a directional MPS reconstruction method that produces results that are more similar to the real core. In 2016, Zuo et al. [172] proposed a method to improve the reconstruction speed of the MPS algorithm by using the half-template technique. In 2017, Peng [173] established a digital core of a shale reservoir by using the MPS method and concluded that the MPS method can effectively reconstruct large-scale digital cores of shale reservoirs by evaluating the function and permeability results. In 2018, Wu et al. [149] combined the MPS method with 3D digital core data to reconstruct a more accurate digital core.Generally, current research on MPS reconstruction methods focuses mainly on three aspects: accelerating the reconstruction speed, reducing the required memory for reconstruction, and using multiscale templates to improve the reconstruction accuracy. Digital cores that are reconstructed with the MPS method display good connectivity; hence, the MPS approach is a widely used reconstruction method. However, similar to the MCMC reconstruction method, we must study and discuss the effectiveness of reconstructed rocks with complicated pore structures and fine-scale descriptions. In addition, the stochastic reconstruction method based on the statistical constraints of 2D slice morphological characteristics does not consider the mathematical statistics of 2D slices; consequently, 3D digital cores that are constructed with this method have some hidden dangers, and certain statistics may not reflect the characteristics of 2D slices. Compared to other stochastic simulation methods, the greatest advantage of MPS is that it can effectively replicate 2D or 3D models of the pore structure and reconstruct the long-distance connectivity of a pore space. Thus, digital cores that are constructed with the MPS technique are regarded as having the best long-distance connectivity. Moreover, this technique can better describe the shape of the pore space, which is difficult to achieve by other reconstruction methods. Unfortunately, similar to other methods based on the reconstruction of morphological features, the MPS approach has a common problem wherein the training image must be representative. In addition, this technique suffers from a relatively slow computational speed. Overall, MPS generates reliable results compared to other reconstruction methods based on the pore morphology. Since its proposal, the MPS digital core-reconstruction method has been employed by a large number of scholars and has continued to improve.Accordingly, the MPS method is currently the most common technique for reconstructing complex 3D digital cores with 2D slices [174]. The reconstructed digital core is more accurate, but the algorithm itself is difficult to solve with the stochastic reconstruction method based on the statistical constraints of 2D slice mathematical features because of the fixed nature of the algorithm. Hence, the stochastic reconstruction method based on the statistical constraints of 2D slice morphological characteristics is mainly used to improve the reconstruction speed and research multiscale reconstruction.


### 3.3. Physics Process-Based Stochastic Reconstruction Method

Stochastic reconstruction methods based on physical processes, that is, methods that consider the sedimentation and diagenesis of real rocks and stochastically simulate these processes, are different from the other two reconstruction methods that were mentioned above. Common examples of such reconstruction methods include random particle stacking considering gravity and random particle stacking based on the discrete element method (DEM).In 1992, Bryant et al. [175] proposed the reconstruction of digital cores with the accumulation of isospheres to reflect geological formation processes; their method was later improved in 1996 [176]. Subsequently, Bakke and Øren [177–179] established 3D digital cores that used spheres of different particle sizes and clearly delineated the deposition, compaction, and diagenesis processes. Through the correlation of two points, the use of local porosity distribution functions, the calculation of the local seepage probability function, and the simulation of formation factors, the reconstructed digital cores matched the digital cores that were obtained through CT scanning. The random particle stacking method can be divided into four steps. First, the particle size-distribution curve can be constructed by considering the number of particles in the 2D slice, the particle size, and the particle size distribution. Second, particles in free fall are randomly selected for deposition simulation in a fixed box; that is, no lateral force is assumed, the particles that fall to a stable position are measured, and the particles fall until the stopping requirement is met. Third, the compaction process is computed by using a degree-of-compaction formula and compaction factors. Finally, the cementing factor is used to simulate the cementing process with a degree-of-cementing formula [180]. Figure 12 shows the geometry process in the physics process-based stochastic reconstruction method, and Figure 13 shows a 3D digital core that was reconstructed with the physics-based method.Figure 12
Geometry process in physics process-based stochastic reconstruction method. (a) Shows 3D digital core results before meshing, and (b) shows 3D digital core results after meshing. Figure from [63].
(a)
(b)Figure 13
The 3D digital sandstone cores were reconstructed using a process-based reconstruction method. Among them, blue is the skeleton part, and the transparent part is the pore part. Pore space connectivity of 3D digital cores is very good. However, it can also be seen that due to the lack of detailed description of the diagenesis of the rock, the three-dimensional digital core obtained by the process method has a very single pore space and is difficult to reconstruct the complicated pore structure.In 2003, Jin et al. suggested that sedimentation should consider the effects of gravity, contact forces, and resistance; accordingly, these authors reconstructed a digital core by using a DEM. This method yields 3D digital cores with a smaller porosity and permeability [181]. In 2012, Zhu et al. [182] proposed a digital core-reconstruction method by changing the shape of an irregular particle. In 2013, Yan et al. [183] compared digital cores that were reconstructed with the physics-based method with artificial cores and concluded that the numerical cores from this method were similar to the artificial cores with regard to the results of the evaluation functions. Zhao et al. [184] conducted a study based on the physics-based process law and achieved good results. In 2015, Zou et al. [185] utilized a DEM to reconstruct 3D digital cores in tight sandstone reservoirs and found these cores to be similar to real digital cores according to NMR simulations. In 2017, Zhu and Yu [186] improved the diagenesis process in this method by converting the voxels adjacent to the framework to build cements. In 2018, Tian et al. [187] proposed a new physics-based process approach that considers the physical influences of gravity and particle collision on a mineral and improves the steady-state particle position search algorithm; thus, deposited particles can quickly converge to stable positions. Their process method improved the disadvantages of previous physics-based process methods, which cannot simulate more complex rocks to a certain extent, and extended its scope of application.Stochastic reconstruction methods based on physical processes have obvious advantages and disadvantages. One advantage is that the obtained rock connectivity is better for simulating real geological processes and is more suitable for simulating the simple diagenesis of rocks than other methods. This advantage is critical, especially for rock seepage and conductivity studies, because differences in the connectivity can have very prominent effects.Among the disadvantages of these methods, rocks such as carbonates and shales that have a complex diagenetic history cannot be reconstructed. As mentioned earlier, process methods can better simulate the rock deposition process, which is similar regardless of whether the rock structure is complex. However, the simulation of diagenesis in the physics-based process method is too simple because it uses only one formula to represent the diagenetic process, which is obviously unreasonable because the diagenetic process should be more complex than the deposition process. The next step in research on process simulation methods should address diagenetic processes, the study of which first must consider the different diagenetic characteristics of different rock types. For example, carbonate rocks are easily eroded [188] and are also prone to fractures and pressure fractures [189, 190], whereas shale reservoirs must consider the conversion of organic matter [191] and the formation of microfractures [192], and corresponding algorithms must be established. A hybrid method can be used to strengthen the simulation of diagenesis, and the targeted simulation of rock diagenesis should improve the use of physics-based process simulation methods for complex rocks.


### 3.4. Hybrid Stochastic Reconstruction Method

According to the above summary of 2D slice-based 3D digital core-reconstruction methods, digital cores that are constructed with different methods evidently produce different results. All these methods have unique advantages and disadvantages, and their merits are difficult to evaluate. A method must be developed to select the technique to be used to reconstruct 3D digital cores. Table3 compares the advantages and disadvantages of common reconstruction methods.Table 3
Common reconstruction methods contrast.
Reconstruction method
Advantage
Disadvantages
Stochastic reconstruction method based on 2D slice statistical constraint features
Gaussian field simulation method
Reconstruction speed is high
Connectivity of reconstructed 3D digital core is poor
Stochastic search algorithm based on majority operator
Good connectivity
Connectivity decreases with decreasing porosity
Simulated annealing method
Allow to join as many evaluation function constraints
The accuracy is limited by the evaluation effect of the evaluation function
Sequence indicator simulation method
The evaluation results for the constraints are similar to the actual core results
Reconstructed 3D digital core connectivity is poor
Stochastic reconstruction method based on 2D slice morphological characteristics statistical constraint
MCMC simulation method
Reconstructed 3D digital cores feature similarities to the original data
Reconstruction results have a certain degree of randomness
MPS simulation method
Good connectivity, pore morphological characteristics similar to the actual rock
Slow calculation, memory footprint
Physics process-based stochastic reconstruction method
Physics process simulation method
Connectivity is very good
Cannot reconstruct complex rocks with diagenesisBased on Table3, different reconstruction algorithms have different emphases in the reconstruction of 3D digital cores. In response to these features, relevant scholars have proposed a hybrid method that combines different reconstruction methods.When reconstructing digital cores, the differences between methods must be understood to build a more accurate digital core. Hidajat et al. [193] proposed a 3D digital core-reconstruction method that uses the Gaussian field method combined with the SA method: the former is used to reconstruct the initial 3D digital core, which is then used to optimize the SA algorithm. The results from their study showed that the 3D digital cores that were obtained with this hybrid method were similar to real 3D digital cores. Liu et al. [194, 195] proposed a 3D digital core-reconstruction method that uses a combination of the physics-based process simulation method and SA method; the latter was used to improve the simulation effect of diagenesis in the reconstructed 3D digital core. The results showed that the connectivity within the 3D digital core from the hybrid reconstruction system was more realistic than the reconstructed 3D digital core from the SA method (Figure 14). The simulation results from this method were the same as the actual experimental results, and the reconstruction speed was faster than the original speed.Figure 14
Compared with the experimental data, the three-dimensional digital core reconstructed by the hybrid method has good correspondence with the experimental results in the high-porosity part. The three-dimensional digital core constructed by this method is more realistic. Figure reprinted from [195].In 2016, Mo et al. [104] combined the SA algorithm with the stochastic search algorithm based on the majority operator to propose a supplementary SA optimization scheme for 3D digital cores. This scheme improves the shape and pore connectivity of the pore spaces in reconstructed digital cores, and the constructed 3D digital cores are more realistic.Nevertheless, most scholars have overlooked the potential of hybrid reconstruction methods. Through the above analysis of current mainstream reconstruction methods, each method evidently has obvious advantages and disadvantages. Consequently, whether one method can completely replace the others is difficult to determine, and no universal digital core-reconstruction method exists. Therefore, further research should be conducted on combining different methods to greatly improve the accuracy of reconstructed digital cores. Regrettably, few studies have focused on this aspect to date.Regarding hybrid methods, a physics-based process method must be combined with another method if high pore connectivity is required. For reservoirs with greater primary porosity, such as medium- to high-porosity and high-permeability sandstone reservoirs, the pore connectivity is usually better, and the degrees of sorting and rounding are relatively consistent. For this purpose, the process method + optimization method is the most suitable combination. For reservoirs with greater secondary porosity, the pore morphology is not as uniform as that of a primary porosity-based reservoir. Clearly, the pore morphology of the reservoir is very important because a pore morphology that consists of a higher abundance of secondary pores, such as pinholes and cracks in carbonates, is related to the mineral composition, secondary effects, etc. In this case, one should use a reconstruction method that considers the pore morphology and then optimize the digital core based on the actual situation to improve the integrity of the reconstructed core. In contrast, a technique that is similar to the majority operator method, which is suitable for fine-tuning a digital core, is more suitable for the final adjustment of a digital core and the final stage of reconstructing a digital rock core. In summary, hybrid random reconstruction methods, process methods, and reconstruction methods based on the pore morphology are suitable for initially modelling digital cores, whereas optimization-based methods are suitable for further adjustments to the initial digital core, and the details can be fine-tuned by using the majority operator method. Thus, hybrid random reconstruction methods are very important.Generally, hybrid methods have achieved some good results regarding the reconstruction of 3D digital cores with different lithologies. The next step is to accurately identify the differences between these methods and create a reasonable combination of methods to both speed up the model reconstruction and improve the effectiveness of model reconstruction. Different combinations of models have some significance for future research.


## 3.1. Stochastic Reconstruction Based on Statistics

### 3.1.1. Stochastic Reconstruction Based on 2D Slice Statistical Constraint Features

The stochastic reconstruction method is the most common method for reconstructing 3D digital cores based on one or more slices of 2D slice data, various mathematical-feature statistics, and functional information, such as porosity data. Similar to facial reconstruction with abstract facial features (e.g., eye size), the stochastic reconstruction method uses the abstract mathematical features of extracted 2D slices to build the most accurate 3D digital cores [86–90].(1) Gaussian Field Method. In 1974, Joshi was the first to use stochastic information that was derived from the statistics of 2D slices to randomly reconstruct rocks [91]. Joshi proposed a reconstruction technique known as the Gaussian field method, which can be roughly divided into three steps. First, a Gaussian field is randomly generated, and the Gaussian field is filled with independent Gaussian variables. Then, the variables are linearly transformed to correlate the variables. The constraints that are used during the transformation are the core porosity values and the two-point correlation function. Finally, the Gaussian field from the second step is transformed into a digital core by using a nonlinear transformation. However, Joshi built only 2D digital core models at the time.In 1984, Quiblier [92] improved upon the primitive Gaussian field-reconstruction method and constructed 3D digital cores. Subsequently, Adler et al. [93] introduced periodic boundary conditions and the Fourier transform to improve the speed and efficiency of the model, respectively. Furthermore, Hilfer [94] added porosity and probability distribution functions to the model to improve its accuracy. Torquato and Lu [95] used the chord-length distribution function to describe pore space features. Ioannidis et al. [96] added the fast Fourier transform to improve the microscopic pore structure accuracy of the constructed digital core. Giona and Adrover [97] used linear filters to improve the accuracy of the reconstructed model. Roberts [98] replaced the two-point probability function with a multipoint probability function. Most recently, Liang et al. [99] improved the Fourier transform by using a truncated Gaussian method.The abovementioned methods clearly improved the use of the Gaussian field reconstruction method; however, this technique has been criticized because of the poor connectivity within the constructed digital core. The connectivity of the reconstructed 3D digital core does not greatly improve with the use of any of the abovementioned methods; thus, further development of the Gaussian field reconstruction method is necessary. Generally, the advantage of the Gaussian field method is its fast modelling speed. If the accuracy of the digital core is not sufficient, then no quantitative calculation is required if the core to be reconstructed is a high-porosity rock. If a more accurate digital core is required, then the number of constraints must be increased; that is, the evaluation function must be enhanced. However, if additional constraints are included, the modelling speed of the Gaussian field method diminishes, so long-distance connectivity problems cannot be solved, and the accuracy becomes insufficient. To obtain a more accurate 3D digital core in this scenario, one should use a digital 3D core reconstruction method based on an optimization algorithm (such as the simulated annealing (SA) algorithm) with similar principles. Accordingly, the Gaussian field reconstruction approach may be replaced eventually.(2) Stochastic Search Algorithm Based on a Majority Operator. In 2007, Zhao et al. proposed a random search method based on a majority operator and evaluated the permeability of the reconstructed digital core [100]. The main modelling process of the majority operator is divided into four steps. First, the 2D slice two-point probability and linear path functions are calculated. Second, datasets that are the same size as the 3D digital cores are randomly generated (for example, a 300×300×300 voxel digital core is randomly generated when reconstructing a 300×300×300 voxel digital core), and the difference between the two-point probability and linear path functions is calculated as the objective function value. Third, the randomly generated 3D digital cores are optimized based on the majority operator-based random search strategy. Seven main models are used for the optimizing operators: D2Q9, D3Q13, D3Q15, D3Q19, D1Q3, D2Q5, and D3Q27 [101]. Assuming that the D3Q27 model is selected, a voxel that points to a hole can be randomly selected in a digital core; then, the binary states of its 26 neighbourhoods (i.e., the porosity or rock skeleton) can be judged. If a point represents a pore among the 26 neighbourhoods that exceeds a certain value (usually set as a constant), the point is selected again; otherwise, the point is selected and replaced with a skeleton point. After each replacement is completed, the target function value is recalculated until the target function value falls below a given value or the target function value no longer changes over several iterations.Only scarce literature is available regarding the multiple-operator random reconstruction method, which is more suitable for the construction of 3D digital cores with large pores (usually more than 20% porosity) [102] because this approach suffers from many weaknesses. First, this method is optimized with a fixed law, which is very slow, and thus easily falls into extreme values, making it difficult for the model to continue the optimization (i.e., SA and other metaheuristics are random and are therefore the best approach to solve the traditional optimization problem). Second, the adjustment of this method is partial, not integral. This partial adjustment can guarantee only the details, while the overall long-range connectivity is difficult to guarantee. Although the majority operator method works better than the Gaussian field method, it still does not address the holistic, random, and long-range connectivity that is required for optimized digital core-reconstruction methods. However, the multiple-operator random reconstruction method can be used as an optimization solution for reconstructing digital cores (described in Section 3.3).(3) Simulated Annealing (SA) Method. In 1997, Hazlett [103] first proposed a new digital core-reconstruction method based on the SA optimization algorithm. After evaluating the 2D slice porosity, the two-point probability function, and the linear path, a 3D digital core that has the same porosity as the 2D slice of the rock is first constructed, after which the digital core is optimized by using an SA optimization algorithm. During each iteration of the optimization, a pore voxel point and skeleton voxel point are randomly selected and exchanged, and the objective function is calculated. If the value of the function decreases, the 3D digital core is updated, the iteration is continued until it terminates, and a reconstructed 3D digital core is obtained. Figure 5 shows the flow chart for reconstructing a 3D digital core with the SA method. Figure 6 shows a 3D digital sandstone core that was obtained with the SA optimization algorithm in Figure 6, in which blue colors denote pores and red colors represent the skeleton.Figure 5
SA method reconstruction of 3D digital core reconstruction.Figure 6
Reconstruction of 3D digital cores by the SA method. The blue part is the pore; the red part is the skeleton. It can be seen that the constructed 3D digital core has better pore connectivity and is similar to the real pore space. Figure reprinted from [104].Subsequently, Yeong, Torquato and Manwart [105–107] proved that the SA method can approximate the evaluation function. Eschricht et al. [108] improved the SA method and the accuracy of the reconstructed model. After building digital cores with the SA method, Zhao and Yao [109] suggested that good evaluation parameters can improve the reconstruction accuracy and consequently proposed an improved method to speed up the model optimization process. Considering the disadvantages of the SA optimization method, Teng et al. [110] proposed a 3D image reconstruction method based on particle swarm optimization (PSO); an improved reconstruction was obtained compared to the Gaussian field reconstruction method and the SA method. Ju et al. [111] proposed an improved SA method that enhanced the modelling effect. Likewise, Zou et al. [112] improved the SA method and proposed a very fast SA optimization method, which was applied to reconstruct shale reservoirs. Frączek et al. [113] proposed a multiscale SA reconstruction method, and Zeng et al. used the SA reconstruction method to reconstruct gas-hydrate digital cores, achieving very good results. Lin et al. [114] proposed a multiscale SA reconstruction method and reconstructed digital carbonate-rock cores with good results. In 2018, Capek proposed an improved SA method to further improve the reconstruction accuracy [115].The effectiveness of reconstructions from the SA method mainly depends on the efficiency of the algorithm and the use of an appropriate evaluation function. With the same evaluation function, the effect of the construction of the SA method is better than that of the Gaussian field method as well as that of the majority operator. The key issue at present is that strong global optimization is the most important component to improve the long-range connectivity regardless of the optimization model. Thus, a large number of evaluation functions that reflect the different details of digital cores can be added to improve the prediction effect. Among the several methods that are currently employed for this task, metaheuristic optimization algorithms, including SA optimization algorithms, are the most promising because of their high optimization efficiency and certain randomness (thereby avoiding extreme values and improving the final optimization effect). In Section4, we also mention that the prediction may be further improved if a new metaheuristic optimization algorithm is used. In addition, the combination of multiple optimization algorithms based on their unique characteristics can further improve the digital core-reconstruction effect.(4) Discussion of the Statistically Constrained Stochastic Reconstruction of Mathematical Features. Many methods are based on the statistically constrained stochastic reconstruction of mathematical features. In this paper, four common methods are introduced. Among them, the most studied is the SA method, which is a reconstruction method that we consider to be very promising. The digital cores that are reconstructed by the SA method boast better connectivity than reconstructed 3D digital cores from the stochastic search algorithm based on the majority operator and Gaussian field methods; thus, the constraint function of the evaluation function can be better executed. In this section, we focus on a discussion of SA methods. To reconstruct a digital core, the effect of an algorithm depends on the required time for the reconstruction and the approximation degree (error size) of the evaluation function. To improve these two components, an appropriate initial model and the algorithm itself should be reasonably improved. The appropriate initial model can be based on a combination of different methods, which are discussed in detail in Section 3.3.The SA algorithm, which was first proposed in 1953, is a classic optimization algorithm that regards the accuracy of the 3D digital core reconstruction as an optimization problem. SA algorithm optimization works very well for rocks with high porosity and high permeability and rocks with uncomplicated pore structures. However, the speed and accuracy are greatly reduced when reconstructing large-scale rocks with complex pore structures [116], indicating that the classical SA algorithm cannot meet complex optimization requirements. Many scholars that have engaged in 3D digital core-reconstruction research have not addressed any improvements to this algorithm itself. Instead, scholars who have focused on optimization research have made improvements to various individual aspects of the SA algorithm, such as the ease with which the algorithm falls into local minima [117, 118], its long execution time [119–125], and the sensitivity of its parameters [126–128]. Since the concept of the metaheuristic optimization algorithm [103] was proposed, research progress regarding this algorithm has been very rapid. Numerous new metaheuristic optimization algorithms based on various phenomena have been proposed. Compared to these metaheuristics, the SA algorithm is primitive. Table 2 summarizes some of the classic metaheuristic optimization algorithms and the latest metaheuristic optimization algorithms. Figure 7 shows a 3D digital core that was reconstructed by PSO and the SA algorithm with the same evaluation function.Table 2
Some of the metaheuristic optimization algorithms.
Imitated object
Representative algorithms
The evolution of nature
Genetic algorithm [129], differential evolution algorithm [130]
Human behaviour
Immune algorithm [131], brainstorm algorithm [132], fireworks algorithm [133]
Physical characteristics
Simulated annealing algorithm [134], intelligent water drop algorithm [135], binary black hole algorithm [136]
Animal behaviour
Particle swarm optimization algorithm [137], ant algorithm [138], bat algorithm [139], krill herd algorithm [140], whale optimization algorithm [141]
Plant behaviour
Invasive weed algorithm [142], rain forest algorithm [143], root growth model [144]Figure 7
Particle swarm optimization (a), Fourier transform-based Gaussian field method (b), and simulated annealing method (c) reconstruction results of 3D digital core. Among them, the 3D digital core with the same background color is the skeleton, and the white part is the pores. It can be observed that the reconstruction effect of particle swarm optimization is better than that of the simulated annealing algorithm, and the obtained 3D digital core has better connectivity in the pore space. Figure from [109].
(a)
(b)
(c)Table2 demonstrates the evident diversity of heuristic optimization algorithms; moreover, only a small number of the many available optimization algorithms are shown in Table 2. Hundreds of optimization algorithms can be used, even if the improved algorithm is not calculated; in addition, no improved algorithms based on the abovementioned metaheuristic optimization algorithm are included in this table. Many of these optimization algorithms have been proven superior to the 35-year-old SA method in terms of the search speed and optimization effect. For the reconstruction of 3D digital cores, one should use a variety of the most recently developed optimization algorithms to determine which is the most suitable for digital core-reconstruction problems and research and to improve the search ability of each algorithm as well as its optimization speed in order to satisfy the needs of higher evaluation functions. This approach is also the best to improve the accuracy of the stochastic reconstruction method based on the statistical constraint of mathematical features. An interdisciplinary method should be used to improve the reliability of a reconstructed 3D digital core. Alternatively, we can combine the characteristics of multiple optimization methods to jointly reconstruct a digital core. For example, the Gaussian field method is used first to optimize the initial stochastic model; at this time, fewer evaluation functions are selected. The results are then transferred to the optimization algorithm, and more evaluation functions are added to refine the digital core modelling process, thereby increasing the long-range connectivity and reliability of the digital core. Finally, a multituning method is used to further fine-tune the model, and then the model more closely resembles the actual rock. This approach can also be regarded as a hybrid optimization method. In summary, the advantage of the 2D slice method based on the statistical constraint of mathematical features is that it considers the integrity of the digital core during its reconstruction; in addition, this technique first generates the 3D random data and then performs the optimization. However, the disadvantages of this approach are obvious. Morphological characteristics such as sex are obviously insufficient and depend on mathematical algorithms.


## 3.1.1. Stochastic Reconstruction Based on 2D Slice Statistical Constraint Features

The stochastic reconstruction method is the most common method for reconstructing 3D digital cores based on one or more slices of 2D slice data, various mathematical-feature statistics, and functional information, such as porosity data. Similar to facial reconstruction with abstract facial features (e.g., eye size), the stochastic reconstruction method uses the abstract mathematical features of extracted 2D slices to build the most accurate 3D digital cores [86–90].(1) Gaussian Field Method. In 1974, Joshi was the first to use stochastic information that was derived from the statistics of 2D slices to randomly reconstruct rocks [91]. Joshi proposed a reconstruction technique known as the Gaussian field method, which can be roughly divided into three steps. First, a Gaussian field is randomly generated, and the Gaussian field is filled with independent Gaussian variables. Then, the variables are linearly transformed to correlate the variables. The constraints that are used during the transformation are the core porosity values and the two-point correlation function. Finally, the Gaussian field from the second step is transformed into a digital core by using a nonlinear transformation. However, Joshi built only 2D digital core models at the time.In 1984, Quiblier [92] improved upon the primitive Gaussian field-reconstruction method and constructed 3D digital cores. Subsequently, Adler et al. [93] introduced periodic boundary conditions and the Fourier transform to improve the speed and efficiency of the model, respectively. Furthermore, Hilfer [94] added porosity and probability distribution functions to the model to improve its accuracy. Torquato and Lu [95] used the chord-length distribution function to describe pore space features. Ioannidis et al. [96] added the fast Fourier transform to improve the microscopic pore structure accuracy of the constructed digital core. Giona and Adrover [97] used linear filters to improve the accuracy of the reconstructed model. Roberts [98] replaced the two-point probability function with a multipoint probability function. Most recently, Liang et al. [99] improved the Fourier transform by using a truncated Gaussian method.The abovementioned methods clearly improved the use of the Gaussian field reconstruction method; however, this technique has been criticized because of the poor connectivity within the constructed digital core. The connectivity of the reconstructed 3D digital core does not greatly improve with the use of any of the abovementioned methods; thus, further development of the Gaussian field reconstruction method is necessary. Generally, the advantage of the Gaussian field method is its fast modelling speed. If the accuracy of the digital core is not sufficient, then no quantitative calculation is required if the core to be reconstructed is a high-porosity rock. If a more accurate digital core is required, then the number of constraints must be increased; that is, the evaluation function must be enhanced. However, if additional constraints are included, the modelling speed of the Gaussian field method diminishes, so long-distance connectivity problems cannot be solved, and the accuracy becomes insufficient. To obtain a more accurate 3D digital core in this scenario, one should use a digital 3D core reconstruction method based on an optimization algorithm (such as the simulated annealing (SA) algorithm) with similar principles. Accordingly, the Gaussian field reconstruction approach may be replaced eventually.(2) Stochastic Search Algorithm Based on a Majority Operator. In 2007, Zhao et al. proposed a random search method based on a majority operator and evaluated the permeability of the reconstructed digital core [100]. The main modelling process of the majority operator is divided into four steps. First, the 2D slice two-point probability and linear path functions are calculated. Second, datasets that are the same size as the 3D digital cores are randomly generated (for example, a 300×300×300 voxel digital core is randomly generated when reconstructing a 300×300×300 voxel digital core), and the difference between the two-point probability and linear path functions is calculated as the objective function value. Third, the randomly generated 3D digital cores are optimized based on the majority operator-based random search strategy. Seven main models are used for the optimizing operators: D2Q9, D3Q13, D3Q15, D3Q19, D1Q3, D2Q5, and D3Q27 [101]. Assuming that the D3Q27 model is selected, a voxel that points to a hole can be randomly selected in a digital core; then, the binary states of its 26 neighbourhoods (i.e., the porosity or rock skeleton) can be judged. If a point represents a pore among the 26 neighbourhoods that exceeds a certain value (usually set as a constant), the point is selected again; otherwise, the point is selected and replaced with a skeleton point. After each replacement is completed, the target function value is recalculated until the target function value falls below a given value or the target function value no longer changes over several iterations.Only scarce literature is available regarding the multiple-operator random reconstruction method, which is more suitable for the construction of 3D digital cores with large pores (usually more than 20% porosity) [102] because this approach suffers from many weaknesses. First, this method is optimized with a fixed law, which is very slow, and thus easily falls into extreme values, making it difficult for the model to continue the optimization (i.e., SA and other metaheuristics are random and are therefore the best approach to solve the traditional optimization problem). Second, the adjustment of this method is partial, not integral. This partial adjustment can guarantee only the details, while the overall long-range connectivity is difficult to guarantee. Although the majority operator method works better than the Gaussian field method, it still does not address the holistic, random, and long-range connectivity that is required for optimized digital core-reconstruction methods. However, the multiple-operator random reconstruction method can be used as an optimization solution for reconstructing digital cores (described in Section 3.3).(3) Simulated Annealing (SA) Method. In 1997, Hazlett [103] first proposed a new digital core-reconstruction method based on the SA optimization algorithm. After evaluating the 2D slice porosity, the two-point probability function, and the linear path, a 3D digital core that has the same porosity as the 2D slice of the rock is first constructed, after which the digital core is optimized by using an SA optimization algorithm. During each iteration of the optimization, a pore voxel point and skeleton voxel point are randomly selected and exchanged, and the objective function is calculated. If the value of the function decreases, the 3D digital core is updated, the iteration is continued until it terminates, and a reconstructed 3D digital core is obtained. Figure 5 shows the flow chart for reconstructing a 3D digital core with the SA method. Figure 6 shows a 3D digital sandstone core that was obtained with the SA optimization algorithm in Figure 6, in which blue colors denote pores and red colors represent the skeleton.Figure 5
SA method reconstruction of 3D digital core reconstruction.Figure 6
Reconstruction of 3D digital cores by the SA method. The blue part is the pore; the red part is the skeleton. It can be seen that the constructed 3D digital core has better pore connectivity and is similar to the real pore space. Figure reprinted from [104].Subsequently, Yeong, Torquato and Manwart [105–107] proved that the SA method can approximate the evaluation function. Eschricht et al. [108] improved the SA method and the accuracy of the reconstructed model. After building digital cores with the SA method, Zhao and Yao [109] suggested that good evaluation parameters can improve the reconstruction accuracy and consequently proposed an improved method to speed up the model optimization process. Considering the disadvantages of the SA optimization method, Teng et al. [110] proposed a 3D image reconstruction method based on particle swarm optimization (PSO); an improved reconstruction was obtained compared to the Gaussian field reconstruction method and the SA method. Ju et al. [111] proposed an improved SA method that enhanced the modelling effect. Likewise, Zou et al. [112] improved the SA method and proposed a very fast SA optimization method, which was applied to reconstruct shale reservoirs. Frączek et al. [113] proposed a multiscale SA reconstruction method, and Zeng et al. used the SA reconstruction method to reconstruct gas-hydrate digital cores, achieving very good results. Lin et al. [114] proposed a multiscale SA reconstruction method and reconstructed digital carbonate-rock cores with good results. In 2018, Capek proposed an improved SA method to further improve the reconstruction accuracy [115].The effectiveness of reconstructions from the SA method mainly depends on the efficiency of the algorithm and the use of an appropriate evaluation function. With the same evaluation function, the effect of the construction of the SA method is better than that of the Gaussian field method as well as that of the majority operator. The key issue at present is that strong global optimization is the most important component to improve the long-range connectivity regardless of the optimization model. Thus, a large number of evaluation functions that reflect the different details of digital cores can be added to improve the prediction effect. Among the several methods that are currently employed for this task, metaheuristic optimization algorithms, including SA optimization algorithms, are the most promising because of their high optimization efficiency and certain randomness (thereby avoiding extreme values and improving the final optimization effect). In Section4, we also mention that the prediction may be further improved if a new metaheuristic optimization algorithm is used. In addition, the combination of multiple optimization algorithms based on their unique characteristics can further improve the digital core-reconstruction effect.(4) Discussion of the Statistically Constrained Stochastic Reconstruction of Mathematical Features. Many methods are based on the statistically constrained stochastic reconstruction of mathematical features. In this paper, four common methods are introduced. Among them, the most studied is the SA method, which is a reconstruction method that we consider to be very promising. The digital cores that are reconstructed by the SA method boast better connectivity than reconstructed 3D digital cores from the stochastic search algorithm based on the majority operator and Gaussian field methods; thus, the constraint function of the evaluation function can be better executed. In this section, we focus on a discussion of SA methods. To reconstruct a digital core, the effect of an algorithm depends on the required time for the reconstruction and the approximation degree (error size) of the evaluation function. To improve these two components, an appropriate initial model and the algorithm itself should be reasonably improved. The appropriate initial model can be based on a combination of different methods, which are discussed in detail in Section 3.3.The SA algorithm, which was first proposed in 1953, is a classic optimization algorithm that regards the accuracy of the 3D digital core reconstruction as an optimization problem. SA algorithm optimization works very well for rocks with high porosity and high permeability and rocks with uncomplicated pore structures. However, the speed and accuracy are greatly reduced when reconstructing large-scale rocks with complex pore structures [116], indicating that the classical SA algorithm cannot meet complex optimization requirements. Many scholars that have engaged in 3D digital core-reconstruction research have not addressed any improvements to this algorithm itself. Instead, scholars who have focused on optimization research have made improvements to various individual aspects of the SA algorithm, such as the ease with which the algorithm falls into local minima [117, 118], its long execution time [119–125], and the sensitivity of its parameters [126–128]. Since the concept of the metaheuristic optimization algorithm [103] was proposed, research progress regarding this algorithm has been very rapid. Numerous new metaheuristic optimization algorithms based on various phenomena have been proposed. Compared to these metaheuristics, the SA algorithm is primitive. Table 2 summarizes some of the classic metaheuristic optimization algorithms and the latest metaheuristic optimization algorithms. Figure 7 shows a 3D digital core that was reconstructed by PSO and the SA algorithm with the same evaluation function.Table 2
Some of the metaheuristic optimization algorithms.
Imitated object
Representative algorithms
The evolution of nature
Genetic algorithm [129], differential evolution algorithm [130]
Human behaviour
Immune algorithm [131], brainstorm algorithm [132], fireworks algorithm [133]
Physical characteristics
Simulated annealing algorithm [134], intelligent water drop algorithm [135], binary black hole algorithm [136]
Animal behaviour
Particle swarm optimization algorithm [137], ant algorithm [138], bat algorithm [139], krill herd algorithm [140], whale optimization algorithm [141]
Plant behaviour
Invasive weed algorithm [142], rain forest algorithm [143], root growth model [144]Figure 7
Particle swarm optimization (a), Fourier transform-based Gaussian field method (b), and simulated annealing method (c) reconstruction results of 3D digital core. Among them, the 3D digital core with the same background color is the skeleton, and the white part is the pores. It can be observed that the reconstruction effect of particle swarm optimization is better than that of the simulated annealing algorithm, and the obtained 3D digital core has better connectivity in the pore space. Figure from [109].
(a)
(b)
(c)Table2 demonstrates the evident diversity of heuristic optimization algorithms; moreover, only a small number of the many available optimization algorithms are shown in Table 2. Hundreds of optimization algorithms can be used, even if the improved algorithm is not calculated; in addition, no improved algorithms based on the abovementioned metaheuristic optimization algorithm are included in this table. Many of these optimization algorithms have been proven superior to the 35-year-old SA method in terms of the search speed and optimization effect. For the reconstruction of 3D digital cores, one should use a variety of the most recently developed optimization algorithms to determine which is the most suitable for digital core-reconstruction problems and research and to improve the search ability of each algorithm as well as its optimization speed in order to satisfy the needs of higher evaluation functions. This approach is also the best to improve the accuracy of the stochastic reconstruction method based on the statistical constraint of mathematical features. An interdisciplinary method should be used to improve the reliability of a reconstructed 3D digital core. Alternatively, we can combine the characteristics of multiple optimization methods to jointly reconstruct a digital core. For example, the Gaussian field method is used first to optimize the initial stochastic model; at this time, fewer evaluation functions are selected. The results are then transferred to the optimization algorithm, and more evaluation functions are added to refine the digital core modelling process, thereby increasing the long-range connectivity and reliability of the digital core. Finally, a multituning method is used to further fine-tune the model, and then the model more closely resembles the actual rock. This approach can also be regarded as a hybrid optimization method. In summary, the advantage of the 2D slice method based on the statistical constraint of mathematical features is that it considers the integrity of the digital core during its reconstruction; in addition, this technique first generates the 3D random data and then performs the optimization. However, the disadvantages of this approach are obvious. Morphological characteristics such as sex are obviously insufficient and depend on mathematical algorithms.


## 3.2. Stochastic Reconstruction Method Based on the Statistical Constraints of 2D Slice Morphological Characteristics

A stochastic reconstruction method based on the statistical constraints of 2D slice morphological characteristics can be used to calculate the morphological features of 2D slices (specifically, the relationships between the positions and relative positions of voxels). This method can also be used for 3D reconstruction purposes based on the morphological features of 2D slices, but it differs from the stochastic method based on the statistical constraints of mathematical features, which considers the correlations between voxels with regard to the connectivity and morphological features of pores. Representative reconstruction algorithms include the Markov chain Monte Carlo (MCMC) reconstruction method and the stochastic reconstruction method based on the statistical constraints of slice morphology characteristics; these algorithms are used to calculate the morphological features of 2D slices (specifically, the relationships between the positions and relative positions of voxels). A reconstructed 3D digital core is then built based on the morphological features of these 2D slices. This approach is different from the stochastic method based on statistically constrained mathematical features, which considers the correlations between voxels with regard to the connectivity and morphological features of pores. Representative reconstruction algorithms include the MCMC reconstruction method and multipoint geostatistics (MPS) [145].


### 3.2.1. Sequence Indicator Simulation Method

In 2003, Keehm [146] proposed a new reconstruction method, namely, a 3D digital core stochastic reconstruction method based on the sequence indicator simulation method from geostatistics. Sequence indicator simulation is a method that indicates the probability distribution field of data by using a directed kriging interpolation technique in combination with a conditional stochastic simulation. In this method, the variogram and porosity values of the 2D slice are first calculated. Then, the path of each network node is randomly assessed, and kriging interpolation is applied to the indicator variable to determine the probability that this variable is a discrete variable at the node. Then, the discrete variable order is determined, and a random number is generated, which determines the random variable type at that point. The nodes are repeatedly updated until each node is simulated, and a 3D digital core is ultimately obtained. In 2007, Zhu and Tao [147] reconstructed a 3D digital core by using cast sheet data. The variation function of the reconstructed digital core is considered to correspond well with the variation function that is calculated from the reconstructed casting slice. The sequence indicator simulation method produces reconstructed 3D digital cores that are very similar to real digital cores, so this method can be applied to the simulation of electrical and elastic rock properties. In 2011, Liu and Mu [148] proposed an improved sequence indicator simulation method and applied it to calculate the permeability of a digital core.Unfortunately, the sequence indicator simulation method does not fundamentally solve the reconstructed digital core connectivity problem. The greatest advantage of the sequence indicator simulation method is that it has certain integrity and correlations between voxels, which is very important for the reconstruction of digital cores to ensure long-range connectivity. The relevance of each voxel is greatly considered. In contrast, the weakness of the sequence indicator simulation method is that the rules to establish the digital core are directly fixed from the image, which causes two hidden problems. First, when using this method for the reconstruction of a digital core, the reconstruction accuracy is directly derived from the image. For images with very high heterogeneity, the image representation is weak, and the reconstruction effect is influenced accordingly. Second, because the evaluation function is relatively fixed, developing it further is difficult, so the model has weaker potential. The above problem is a common problem of modelling methods that are based on the morphological features of images, but the disadvantages of the sequence indicator simulation method are more serious. Compared to the MCMC and MPS methods, the digital cores that are reconstructed by the sequence indicator simulation method are also less reliable and are especially weaker overall than those that are reconstructed by the MPS method, which considers several consecutive voxel points in the vicinity and thereby more fully considers the digital core. The subsequent implementation of the simulation method is significantly worse than the implementation of both the MCMC and MPS methods. In summary, the sequence indicator simulation method does not fundamentally solve the reconstructed digital core connectivity problem.


### 3.2.2. Markov Chain Monte Carlo (MCMC) Method

In 2004, the MCMC reconstruction method was proposed by Wu et al., who applied this technique to the 2D reconstruction of soil structure [149]. In 2006, Wu et al. extended the 2D digital core-reconstruction method into three dimensions and proposed an MCMC-based 3D digital core-reconstruction method [150]. This technique introduces the Markov chain and traverses all two-point and five-point neighbourhoods in a 2D slice image to calculate the conditional probability of the neighbourhood template; this method also determines the state of each point in the reconstructed image based on conditional probabilities. Markov chains describe a sequence of states (analogous to the notion of a template). The state value of each position in the sequence depends on the state of the previous position, and the probability of describing this state is called the transition probability [151, 152]. This reconstruction method can be roughly divided into four steps. First, we establish a Markov chain with a stable probability distribution function by means of traversal scanning. When the chain is performed over a sufficient distance in the image, an important statistical feature of the original image is obtained. After the conditional probability is obtained, the image reconstruction can be performed. Then, the first row of the first layer is reconstructed, and the average porosity is determined by using 2D slices in the x, y, and z directions (slices can be reconstructed in three directions from a single 2D slice). From the second voxel, a two-point neighbourhood template is used to reconstruct the three-point neighbourhood system from the third point to reconstruct the first row of voxel values. The third step is to reconstruct the first layer of the 3D digital core. The first point is reconstructed with two-point neighbourhood templates, and the second point is reconstructed by four-point neighbourhood templates. The remaining points are reconstructed by using the five- and six-point neighbourhood systems. Subsequently, the remaining layers of the 3D digital core are reconstructed, and the first layer is reconstructed vertically by using the previous steps. Then, the first voxel of the second layer utilizes seven- and eight-point neighbourhood systems for reconstruction. Finally, the second line of the remaining components of the 3D digital core can be reconstructed with 10- and 11-point neighbourhood systems. A diagram of the MCMC reconstruction method is shown in Figure 8, and Figure 9 displays a digital core of the Berea sandstone that was reconstructed with the MCMC method.Figure 8
Reconstruction process of 3D digital cores using MCMC reconstruction method. (a) Shows reconstruction the first voxel, (b) shows reconstruction the first row, (c) shows reconstruction the first level, and (d) shows reconstruction the remaining layers of the 3D digital core. Figure from [153].
(a)
(b)
(c)
(d)Figure 9
Reconstruction of 3D digital cores using MCMC reconstruction method. Among them, black is the pore part, and white is the skeleton part. The 3D digital cores obtained by the MCMC method have good connectivity of pore space and more concentrated pore space.Wang et al. [154] used the MCMC method to reconstruct the pores of shale rocks in two dimensions. Nie et al. [154] modelled a shale-gas reservoir with a large-scale 3D digital core by using the 3D MCMC method. Guo et al. [155] used a random walk algorithm to simulate the NMR response of a sandstone reservoir in a reconstructed 3D digital core based on the MCMC algorithm. Overall, the MCMC-reconstructed sandstone digital cores match the actual cores, and the NMR response values are similar.The MCMC method is a very reliable digital core-reconstruction technique. The constructed digital core boasts good connectivity. However, MCMC-reconstructed digital cores have weak anisotropy, and the pore throat radius distribution is very concentrated. These factors are the major disadvantages for reservoirs with strong heterogeneity. Compared to the sequential indicator simulation method, the MCMC method considers the voxel configuration around each voxel to be determined and thus obtains a digital core with better connectivity, especially with regard to the distribution of pores, when the core is more uniform. However, the most fundamental cause of the weak anisotropy in reconstructed digital cores is that this hypothesis does not conform to a heterogeneous rock. The assumption of this technique is that the state of any point depends only on the states of a few neighbouring points. Then, after counting the probability function of the training image, the voxels that actually must be determined are determined only by points that are close to each other, whereas the distances to voxel points that are far from the points to be determined are ignored, resulting in the absence of long-distance differences and eventually producing weak anisotropy. Nevertheless, the MCMC method is generally more accurate than the sequential indicator simulation method; consequently, more studies have been performed based on the reconstruction of digital cores with the MCMC method. In conclusion, 3D digital cores that are reconstructed by the MCMC algorithm match real digital cores, but whether these cores are applicable to rocks with very complicated pore structures remains uncertain.


### 3.2.3. Multipoint Geostatistics (MPS) Method

In 2004, Okabe and Blunt proposed a digital core-reconstruction method based on the statistical method of MPS [42]. The MPS method extracts multipoint statistical information from training images (i.e., 2D slices) and reconstructs the template and reconstruction mode; then, this mode is used to perform a 3D digital core reconstruction. Compared to the MCMC method, this method is more biased toward the reconstruction of 3D digital cores with the morphological features of 2D slices. The reconstruction of a digital core with the MPS method is divided into four steps. First, a 2D slice search template and a search tree are established. Second, the raw data are reloaded to the nearest simulation grid node and fixed during the simulation. Next, we define a path that randomly accesses all voxels. We use the search template to define conditional data events, compute the conditional probability distribution function (CPDF) based on the search tree, calculate the simulated values with the function, and generate new 2D images through iteration. Finally, we employ the reconstructed image as a new training image to generate the next layer of the image, which is then used to generate a new 3D digital core. Figure 10 shows a reconfigured flowchart of the MPS reconstruction method, and Figure 11 shows a 3D digital shale-rock core that was reconstructed with the MPS method.Figure 10
MPS reconstruction method flow chart.Figure 11
Reconstruction result of the 3D digital core using MPS method. Among them, the blue part is the pore, and the transparent part is the skeleton. It can be seen that even for shale rock, the long-range connectivity of the pore space of the 3D digital core obtained by the MPS method is still good. This shows that the MPS method can reflect the training image well and obtain the 3D digital core with long-range connectivity. Figure reprinted from [155].In 2004, Daı̈an et al. [156] used MPS to reconstruct digital cores and calculate their porosity and permeability, which matched the actual experimental results. Because the MPS method suffers from a slow reconstruction speed and occupies more memory, Wu et al. [157] proposed a solution to replace this pixel-wise distance calculation with a filter score comparison, improving the reconstruction speed by a factor of 10. In 2009, Wang et al. [158] reconstructed 3D digital cores by using 2D CT images and tested the permeability of the reconstructed 3D digital cores by using the lattice Boltzmann method. Teng et al. [159] proposed that higher multipoint statistical stationarity and lower-scale stationarity are likely to produce results that more closely resemble the real 3D structure. In 2012, Tahmasebi et al. [160] proposed a multipoint statistical method based on cross-correlation function geostatistical methods (CCSIM) and conducted continuous research to improve the method [161–165]. Straubhaar et al. [166] proposed a parallel multipoint geostatistical method based on an improved search tree structure, effectively reducing the memory footprint and increasing the computational speed of the method. In 2013, Du and Zhang [167] used multigrid templates to improve the reconstruction effect and achieved good results. In 2014, Zhang et al. [168] proposed an MPS method based on GPU acceleration. In the same year, Gao et al. [169] proposed an MPS method that uses three-step sampling. Zhang et al. [170] proposed an Isomap-based MPS reconstruction method that uses the Isomap method to reduce the dimensions and redundancy. Liu et al. [171] proposed a directional MPS reconstruction method that produces results that are more similar to the real core. In 2016, Zuo et al. [172] proposed a method to improve the reconstruction speed of the MPS algorithm by using the half-template technique. In 2017, Peng [173] established a digital core of a shale reservoir by using the MPS method and concluded that the MPS method can effectively reconstruct large-scale digital cores of shale reservoirs by evaluating the function and permeability results. In 2018, Wu et al. [149] combined the MPS method with 3D digital core data to reconstruct a more accurate digital core.Generally, current research on MPS reconstruction methods focuses mainly on three aspects: accelerating the reconstruction speed, reducing the required memory for reconstruction, and using multiscale templates to improve the reconstruction accuracy. Digital cores that are reconstructed with the MPS method display good connectivity; hence, the MPS approach is a widely used reconstruction method. However, similar to the MCMC reconstruction method, we must study and discuss the effectiveness of reconstructed rocks with complicated pore structures and fine-scale descriptions. In addition, the stochastic reconstruction method based on the statistical constraints of 2D slice morphological characteristics does not consider the mathematical statistics of 2D slices; consequently, 3D digital cores that are constructed with this method have some hidden dangers, and certain statistics may not reflect the characteristics of 2D slices. Compared to other stochastic simulation methods, the greatest advantage of MPS is that it can effectively replicate 2D or 3D models of the pore structure and reconstruct the long-distance connectivity of a pore space. Thus, digital cores that are constructed with the MPS technique are regarded as having the best long-distance connectivity. Moreover, this technique can better describe the shape of the pore space, which is difficult to achieve by other reconstruction methods. Unfortunately, similar to other methods based on the reconstruction of morphological features, the MPS approach has a common problem wherein the training image must be representative. In addition, this technique suffers from a relatively slow computational speed. Overall, MPS generates reliable results compared to other reconstruction methods based on the pore morphology. Since its proposal, the MPS digital core-reconstruction method has been employed by a large number of scholars and has continued to improve.Accordingly, the MPS method is currently the most common technique for reconstructing complex 3D digital cores with 2D slices [174]. The reconstructed digital core is more accurate, but the algorithm itself is difficult to solve with the stochastic reconstruction method based on the statistical constraints of 2D slice mathematical features because of the fixed nature of the algorithm. Hence, the stochastic reconstruction method based on the statistical constraints of 2D slice morphological characteristics is mainly used to improve the reconstruction speed and research multiscale reconstruction.


## 3.2.1. Sequence Indicator Simulation Method

In 2003, Keehm [146] proposed a new reconstruction method, namely, a 3D digital core stochastic reconstruction method based on the sequence indicator simulation method from geostatistics. Sequence indicator simulation is a method that indicates the probability distribution field of data by using a directed kriging interpolation technique in combination with a conditional stochastic simulation. In this method, the variogram and porosity values of the 2D slice are first calculated. Then, the path of each network node is randomly assessed, and kriging interpolation is applied to the indicator variable to determine the probability that this variable is a discrete variable at the node. Then, the discrete variable order is determined, and a random number is generated, which determines the random variable type at that point. The nodes are repeatedly updated until each node is simulated, and a 3D digital core is ultimately obtained. In 2007, Zhu and Tao [147] reconstructed a 3D digital core by using cast sheet data. The variation function of the reconstructed digital core is considered to correspond well with the variation function that is calculated from the reconstructed casting slice. The sequence indicator simulation method produces reconstructed 3D digital cores that are very similar to real digital cores, so this method can be applied to the simulation of electrical and elastic rock properties. In 2011, Liu and Mu [148] proposed an improved sequence indicator simulation method and applied it to calculate the permeability of a digital core.Unfortunately, the sequence indicator simulation method does not fundamentally solve the reconstructed digital core connectivity problem. The greatest advantage of the sequence indicator simulation method is that it has certain integrity and correlations between voxels, which is very important for the reconstruction of digital cores to ensure long-range connectivity. The relevance of each voxel is greatly considered. In contrast, the weakness of the sequence indicator simulation method is that the rules to establish the digital core are directly fixed from the image, which causes two hidden problems. First, when using this method for the reconstruction of a digital core, the reconstruction accuracy is directly derived from the image. For images with very high heterogeneity, the image representation is weak, and the reconstruction effect is influenced accordingly. Second, because the evaluation function is relatively fixed, developing it further is difficult, so the model has weaker potential. The above problem is a common problem of modelling methods that are based on the morphological features of images, but the disadvantages of the sequence indicator simulation method are more serious. Compared to the MCMC and MPS methods, the digital cores that are reconstructed by the sequence indicator simulation method are also less reliable and are especially weaker overall than those that are reconstructed by the MPS method, which considers several consecutive voxel points in the vicinity and thereby more fully considers the digital core. The subsequent implementation of the simulation method is significantly worse than the implementation of both the MCMC and MPS methods. In summary, the sequence indicator simulation method does not fundamentally solve the reconstructed digital core connectivity problem.


## 3.2.2. Markov Chain Monte Carlo (MCMC) Method

In 2004, the MCMC reconstruction method was proposed by Wu et al., who applied this technique to the 2D reconstruction of soil structure [149]. In 2006, Wu et al. extended the 2D digital core-reconstruction method into three dimensions and proposed an MCMC-based 3D digital core-reconstruction method [150]. This technique introduces the Markov chain and traverses all two-point and five-point neighbourhoods in a 2D slice image to calculate the conditional probability of the neighbourhood template; this method also determines the state of each point in the reconstructed image based on conditional probabilities. Markov chains describe a sequence of states (analogous to the notion of a template). The state value of each position in the sequence depends on the state of the previous position, and the probability of describing this state is called the transition probability [151, 152]. This reconstruction method can be roughly divided into four steps. First, we establish a Markov chain with a stable probability distribution function by means of traversal scanning. When the chain is performed over a sufficient distance in the image, an important statistical feature of the original image is obtained. After the conditional probability is obtained, the image reconstruction can be performed. Then, the first row of the first layer is reconstructed, and the average porosity is determined by using 2D slices in the x, y, and z directions (slices can be reconstructed in three directions from a single 2D slice). From the second voxel, a two-point neighbourhood template is used to reconstruct the three-point neighbourhood system from the third point to reconstruct the first row of voxel values. The third step is to reconstruct the first layer of the 3D digital core. The first point is reconstructed with two-point neighbourhood templates, and the second point is reconstructed by four-point neighbourhood templates. The remaining points are reconstructed by using the five- and six-point neighbourhood systems. Subsequently, the remaining layers of the 3D digital core are reconstructed, and the first layer is reconstructed vertically by using the previous steps. Then, the first voxel of the second layer utilizes seven- and eight-point neighbourhood systems for reconstruction. Finally, the second line of the remaining components of the 3D digital core can be reconstructed with 10- and 11-point neighbourhood systems. A diagram of the MCMC reconstruction method is shown in Figure 8, and Figure 9 displays a digital core of the Berea sandstone that was reconstructed with the MCMC method.Figure 8
Reconstruction process of 3D digital cores using MCMC reconstruction method. (a) Shows reconstruction the first voxel, (b) shows reconstruction the first row, (c) shows reconstruction the first level, and (d) shows reconstruction the remaining layers of the 3D digital core. Figure from [153].
(a)
(b)
(c)
(d)Figure 9
Reconstruction of 3D digital cores using MCMC reconstruction method. Among them, black is the pore part, and white is the skeleton part. The 3D digital cores obtained by the MCMC method have good connectivity of pore space and more concentrated pore space.Wang et al. [154] used the MCMC method to reconstruct the pores of shale rocks in two dimensions. Nie et al. [154] modelled a shale-gas reservoir with a large-scale 3D digital core by using the 3D MCMC method. Guo et al. [155] used a random walk algorithm to simulate the NMR response of a sandstone reservoir in a reconstructed 3D digital core based on the MCMC algorithm. Overall, the MCMC-reconstructed sandstone digital cores match the actual cores, and the NMR response values are similar.The MCMC method is a very reliable digital core-reconstruction technique. The constructed digital core boasts good connectivity. However, MCMC-reconstructed digital cores have weak anisotropy, and the pore throat radius distribution is very concentrated. These factors are the major disadvantages for reservoirs with strong heterogeneity. Compared to the sequential indicator simulation method, the MCMC method considers the voxel configuration around each voxel to be determined and thus obtains a digital core with better connectivity, especially with regard to the distribution of pores, when the core is more uniform. However, the most fundamental cause of the weak anisotropy in reconstructed digital cores is that this hypothesis does not conform to a heterogeneous rock. The assumption of this technique is that the state of any point depends only on the states of a few neighbouring points. Then, after counting the probability function of the training image, the voxels that actually must be determined are determined only by points that are close to each other, whereas the distances to voxel points that are far from the points to be determined are ignored, resulting in the absence of long-distance differences and eventually producing weak anisotropy. Nevertheless, the MCMC method is generally more accurate than the sequential indicator simulation method; consequently, more studies have been performed based on the reconstruction of digital cores with the MCMC method. In conclusion, 3D digital cores that are reconstructed by the MCMC algorithm match real digital cores, but whether these cores are applicable to rocks with very complicated pore structures remains uncertain.


## 3.2.3. Multipoint Geostatistics (MPS) Method

In 2004, Okabe and Blunt proposed a digital core-reconstruction method based on the statistical method of MPS [42]. The MPS method extracts multipoint statistical information from training images (i.e., 2D slices) and reconstructs the template and reconstruction mode; then, this mode is used to perform a 3D digital core reconstruction. Compared to the MCMC method, this method is more biased toward the reconstruction of 3D digital cores with the morphological features of 2D slices. The reconstruction of a digital core with the MPS method is divided into four steps. First, a 2D slice search template and a search tree are established. Second, the raw data are reloaded to the nearest simulation grid node and fixed during the simulation. Next, we define a path that randomly accesses all voxels. We use the search template to define conditional data events, compute the conditional probability distribution function (CPDF) based on the search tree, calculate the simulated values with the function, and generate new 2D images through iteration. Finally, we employ the reconstructed image as a new training image to generate the next layer of the image, which is then used to generate a new 3D digital core. Figure 10 shows a reconfigured flowchart of the MPS reconstruction method, and Figure 11 shows a 3D digital shale-rock core that was reconstructed with the MPS method.Figure 10
MPS reconstruction method flow chart.Figure 11
Reconstruction result of the 3D digital core using MPS method. Among them, the blue part is the pore, and the transparent part is the skeleton. It can be seen that even for shale rock, the long-range connectivity of the pore space of the 3D digital core obtained by the MPS method is still good. This shows that the MPS method can reflect the training image well and obtain the 3D digital core with long-range connectivity. Figure reprinted from [155].In 2004, Daı̈an et al. [156] used MPS to reconstruct digital cores and calculate their porosity and permeability, which matched the actual experimental results. Because the MPS method suffers from a slow reconstruction speed and occupies more memory, Wu et al. [157] proposed a solution to replace this pixel-wise distance calculation with a filter score comparison, improving the reconstruction speed by a factor of 10. In 2009, Wang et al. [158] reconstructed 3D digital cores by using 2D CT images and tested the permeability of the reconstructed 3D digital cores by using the lattice Boltzmann method. Teng et al. [159] proposed that higher multipoint statistical stationarity and lower-scale stationarity are likely to produce results that more closely resemble the real 3D structure. In 2012, Tahmasebi et al. [160] proposed a multipoint statistical method based on cross-correlation function geostatistical methods (CCSIM) and conducted continuous research to improve the method [161–165]. Straubhaar et al. [166] proposed a parallel multipoint geostatistical method based on an improved search tree structure, effectively reducing the memory footprint and increasing the computational speed of the method. In 2013, Du and Zhang [167] used multigrid templates to improve the reconstruction effect and achieved good results. In 2014, Zhang et al. [168] proposed an MPS method based on GPU acceleration. In the same year, Gao et al. [169] proposed an MPS method that uses three-step sampling. Zhang et al. [170] proposed an Isomap-based MPS reconstruction method that uses the Isomap method to reduce the dimensions and redundancy. Liu et al. [171] proposed a directional MPS reconstruction method that produces results that are more similar to the real core. In 2016, Zuo et al. [172] proposed a method to improve the reconstruction speed of the MPS algorithm by using the half-template technique. In 2017, Peng [173] established a digital core of a shale reservoir by using the MPS method and concluded that the MPS method can effectively reconstruct large-scale digital cores of shale reservoirs by evaluating the function and permeability results. In 2018, Wu et al. [149] combined the MPS method with 3D digital core data to reconstruct a more accurate digital core.Generally, current research on MPS reconstruction methods focuses mainly on three aspects: accelerating the reconstruction speed, reducing the required memory for reconstruction, and using multiscale templates to improve the reconstruction accuracy. Digital cores that are reconstructed with the MPS method display good connectivity; hence, the MPS approach is a widely used reconstruction method. However, similar to the MCMC reconstruction method, we must study and discuss the effectiveness of reconstructed rocks with complicated pore structures and fine-scale descriptions. In addition, the stochastic reconstruction method based on the statistical constraints of 2D slice morphological characteristics does not consider the mathematical statistics of 2D slices; consequently, 3D digital cores that are constructed with this method have some hidden dangers, and certain statistics may not reflect the characteristics of 2D slices. Compared to other stochastic simulation methods, the greatest advantage of MPS is that it can effectively replicate 2D or 3D models of the pore structure and reconstruct the long-distance connectivity of a pore space. Thus, digital cores that are constructed with the MPS technique are regarded as having the best long-distance connectivity. Moreover, this technique can better describe the shape of the pore space, which is difficult to achieve by other reconstruction methods. Unfortunately, similar to other methods based on the reconstruction of morphological features, the MPS approach has a common problem wherein the training image must be representative. In addition, this technique suffers from a relatively slow computational speed. Overall, MPS generates reliable results compared to other reconstruction methods based on the pore morphology. Since its proposal, the MPS digital core-reconstruction method has been employed by a large number of scholars and has continued to improve.Accordingly, the MPS method is currently the most common technique for reconstructing complex 3D digital cores with 2D slices [174]. The reconstructed digital core is more accurate, but the algorithm itself is difficult to solve with the stochastic reconstruction method based on the statistical constraints of 2D slice mathematical features because of the fixed nature of the algorithm. Hence, the stochastic reconstruction method based on the statistical constraints of 2D slice morphological characteristics is mainly used to improve the reconstruction speed and research multiscale reconstruction.


## 3.3. Physics Process-Based Stochastic Reconstruction Method

Stochastic reconstruction methods based on physical processes, that is, methods that consider the sedimentation and diagenesis of real rocks and stochastically simulate these processes, are different from the other two reconstruction methods that were mentioned above. Common examples of such reconstruction methods include random particle stacking considering gravity and random particle stacking based on the discrete element method (DEM).In 1992, Bryant et al. [175] proposed the reconstruction of digital cores with the accumulation of isospheres to reflect geological formation processes; their method was later improved in 1996 [176]. Subsequently, Bakke and Øren [177–179] established 3D digital cores that used spheres of different particle sizes and clearly delineated the deposition, compaction, and diagenesis processes. Through the correlation of two points, the use of local porosity distribution functions, the calculation of the local seepage probability function, and the simulation of formation factors, the reconstructed digital cores matched the digital cores that were obtained through CT scanning. The random particle stacking method can be divided into four steps. First, the particle size-distribution curve can be constructed by considering the number of particles in the 2D slice, the particle size, and the particle size distribution. Second, particles in free fall are randomly selected for deposition simulation in a fixed box; that is, no lateral force is assumed, the particles that fall to a stable position are measured, and the particles fall until the stopping requirement is met. Third, the compaction process is computed by using a degree-of-compaction formula and compaction factors. Finally, the cementing factor is used to simulate the cementing process with a degree-of-cementing formula [180]. Figure 12 shows the geometry process in the physics process-based stochastic reconstruction method, and Figure 13 shows a 3D digital core that was reconstructed with the physics-based method.Figure 12
Geometry process in physics process-based stochastic reconstruction method. (a) Shows 3D digital core results before meshing, and (b) shows 3D digital core results after meshing. Figure from [63].
(a)
(b)Figure 13
The 3D digital sandstone cores were reconstructed using a process-based reconstruction method. Among them, blue is the skeleton part, and the transparent part is the pore part. Pore space connectivity of 3D digital cores is very good. However, it can also be seen that due to the lack of detailed description of the diagenesis of the rock, the three-dimensional digital core obtained by the process method has a very single pore space and is difficult to reconstruct the complicated pore structure.In 2003, Jin et al. suggested that sedimentation should consider the effects of gravity, contact forces, and resistance; accordingly, these authors reconstructed a digital core by using a DEM. This method yields 3D digital cores with a smaller porosity and permeability [181]. In 2012, Zhu et al. [182] proposed a digital core-reconstruction method by changing the shape of an irregular particle. In 2013, Yan et al. [183] compared digital cores that were reconstructed with the physics-based method with artificial cores and concluded that the numerical cores from this method were similar to the artificial cores with regard to the results of the evaluation functions. Zhao et al. [184] conducted a study based on the physics-based process law and achieved good results. In 2015, Zou et al. [185] utilized a DEM to reconstruct 3D digital cores in tight sandstone reservoirs and found these cores to be similar to real digital cores according to NMR simulations. In 2017, Zhu and Yu [186] improved the diagenesis process in this method by converting the voxels adjacent to the framework to build cements. In 2018, Tian et al. [187] proposed a new physics-based process approach that considers the physical influences of gravity and particle collision on a mineral and improves the steady-state particle position search algorithm; thus, deposited particles can quickly converge to stable positions. Their process method improved the disadvantages of previous physics-based process methods, which cannot simulate more complex rocks to a certain extent, and extended its scope of application.Stochastic reconstruction methods based on physical processes have obvious advantages and disadvantages. One advantage is that the obtained rock connectivity is better for simulating real geological processes and is more suitable for simulating the simple diagenesis of rocks than other methods. This advantage is critical, especially for rock seepage and conductivity studies, because differences in the connectivity can have very prominent effects.Among the disadvantages of these methods, rocks such as carbonates and shales that have a complex diagenetic history cannot be reconstructed. As mentioned earlier, process methods can better simulate the rock deposition process, which is similar regardless of whether the rock structure is complex. However, the simulation of diagenesis in the physics-based process method is too simple because it uses only one formula to represent the diagenetic process, which is obviously unreasonable because the diagenetic process should be more complex than the deposition process. The next step in research on process simulation methods should address diagenetic processes, the study of which first must consider the different diagenetic characteristics of different rock types. For example, carbonate rocks are easily eroded [188] and are also prone to fractures and pressure fractures [189, 190], whereas shale reservoirs must consider the conversion of organic matter [191] and the formation of microfractures [192], and corresponding algorithms must be established. A hybrid method can be used to strengthen the simulation of diagenesis, and the targeted simulation of rock diagenesis should improve the use of physics-based process simulation methods for complex rocks.


## 3.4. Hybrid Stochastic Reconstruction Method

According to the above summary of 2D slice-based 3D digital core-reconstruction methods, digital cores that are constructed with different methods evidently produce different results. All these methods have unique advantages and disadvantages, and their merits are difficult to evaluate. A method must be developed to select the technique to be used to reconstruct 3D digital cores. Table3 compares the advantages and disadvantages of common reconstruction methods.Table 3
Common reconstruction methods contrast.
Reconstruction method
Advantage
Disadvantages
Stochastic reconstruction method based on 2D slice statistical constraint features
Gaussian field simulation method
Reconstruction speed is high
Connectivity of reconstructed 3D digital core is poor
Stochastic search algorithm based on majority operator
Good connectivity
Connectivity decreases with decreasing porosity
Simulated annealing method
Allow to join as many evaluation function constraints
The accuracy is limited by the evaluation effect of the evaluation function
Sequence indicator simulation method
The evaluation results for the constraints are similar to the actual core results
Reconstructed 3D digital core connectivity is poor
Stochastic reconstruction method based on 2D slice morphological characteristics statistical constraint
MCMC simulation method
Reconstructed 3D digital cores feature similarities to the original data
Reconstruction results have a certain degree of randomness
MPS simulation method
Good connectivity, pore morphological characteristics similar to the actual rock
Slow calculation, memory footprint
Physics process-based stochastic reconstruction method
Physics process simulation method
Connectivity is very good
Cannot reconstruct complex rocks with diagenesisBased on Table3, different reconstruction algorithms have different emphases in the reconstruction of 3D digital cores. In response to these features, relevant scholars have proposed a hybrid method that combines different reconstruction methods.When reconstructing digital cores, the differences between methods must be understood to build a more accurate digital core. Hidajat et al. [193] proposed a 3D digital core-reconstruction method that uses the Gaussian field method combined with the SA method: the former is used to reconstruct the initial 3D digital core, which is then used to optimize the SA algorithm. The results from their study showed that the 3D digital cores that were obtained with this hybrid method were similar to real 3D digital cores. Liu et al. [194, 195] proposed a 3D digital core-reconstruction method that uses a combination of the physics-based process simulation method and SA method; the latter was used to improve the simulation effect of diagenesis in the reconstructed 3D digital core. The results showed that the connectivity within the 3D digital core from the hybrid reconstruction system was more realistic than the reconstructed 3D digital core from the SA method (Figure 14). The simulation results from this method were the same as the actual experimental results, and the reconstruction speed was faster than the original speed.Figure 14
Compared with the experimental data, the three-dimensional digital core reconstructed by the hybrid method has good correspondence with the experimental results in the high-porosity part. The three-dimensional digital core constructed by this method is more realistic. Figure reprinted from [195].In 2016, Mo et al. [104] combined the SA algorithm with the stochastic search algorithm based on the majority operator to propose a supplementary SA optimization scheme for 3D digital cores. This scheme improves the shape and pore connectivity of the pore spaces in reconstructed digital cores, and the constructed 3D digital cores are more realistic.Nevertheless, most scholars have overlooked the potential of hybrid reconstruction methods. Through the above analysis of current mainstream reconstruction methods, each method evidently has obvious advantages and disadvantages. Consequently, whether one method can completely replace the others is difficult to determine, and no universal digital core-reconstruction method exists. Therefore, further research should be conducted on combining different methods to greatly improve the accuracy of reconstructed digital cores. Regrettably, few studies have focused on this aspect to date.Regarding hybrid methods, a physics-based process method must be combined with another method if high pore connectivity is required. For reservoirs with greater primary porosity, such as medium- to high-porosity and high-permeability sandstone reservoirs, the pore connectivity is usually better, and the degrees of sorting and rounding are relatively consistent. For this purpose, the process method + optimization method is the most suitable combination. For reservoirs with greater secondary porosity, the pore morphology is not as uniform as that of a primary porosity-based reservoir. Clearly, the pore morphology of the reservoir is very important because a pore morphology that consists of a higher abundance of secondary pores, such as pinholes and cracks in carbonates, is related to the mineral composition, secondary effects, etc. In this case, one should use a reconstruction method that considers the pore morphology and then optimize the digital core based on the actual situation to improve the integrity of the reconstructed core. In contrast, a technique that is similar to the majority operator method, which is suitable for fine-tuning a digital core, is more suitable for the final adjustment of a digital core and the final stage of reconstructing a digital rock core. In summary, hybrid random reconstruction methods, process methods, and reconstruction methods based on the pore morphology are suitable for initially modelling digital cores, whereas optimization-based methods are suitable for further adjustments to the initial digital core, and the details can be fine-tuned by using the majority operator method. Thus, hybrid random reconstruction methods are very important.Generally, hybrid methods have achieved some good results regarding the reconstruction of 3D digital cores with different lithologies. The next step is to accurately identify the differences between these methods and create a reasonable combination of methods to both speed up the model reconstruction and improve the effectiveness of model reconstruction. Different combinations of models have some significance for future research.


## 4. Discussion of Trends and Challenges

As a new discipline direction, the reconstruction and simulation of a digital core can resolve the problems that are associated with experimental costs and practical difficulties when encountering some tight rocks. This method can aid in the exploration and development of oilfields, which is of great significance. Obtaining accurate 3D digital cores constitutes the premise of accurate simulation experiments and represents the most important research direction regarding digital cores. However, the current modelling methods are far from perfect. Generally, the actual problems that are considered are not meticulous. In this section, we provide a detailed and reasonable view of digital core-reconstruction techniques according to their development and discuss the urgent problems that must be resolved in exploration and development.


### 4.1. Challenges of Reconstructing 3D Digital Cores Based on Physical Experimental Methods

Existing physical experimental methods to reconstruct 3D digital cores encounter fewer challenges than the creation of 3D digital cores based on 2D slices because the former produce more realistic results. Presently, three key problems must be further studied to establish a 3D digital-core method based on physical experimentation. First, the selected rock sample size and optimal scan resolution is important. For a rock with a small pore size, high-resolution instruments can be used to accurately identify the pore space; conversely, a rock with a larger pore size is recommended for low-resolution pore space identification approaches. Second, binary images require further study. The porosity cannot be determined with a single threshold because of obscure skeleton-porosity boundaries in greyscale images. Additionally, the porosity cannot be determined based on physical property experiments, and the microporosity may be difficult to measure through porosity experiments [196–198]. Therefore, the determination of the segmentation threshold is another problem to be addressed. Finally, the size of the 3D digital core must be modest because of computational limitations. Thus, the determination of an appropriate 3D digital-core size for simulation experiments is also meaningful, and the task of building 3D digital cores based on physical experimentation is far less challenging than the reconstruction of 3D digital cores based on 2D slices.


### 4.2. Trends and Challenges of Reconstructing 3D Digital Cores Based on 2D Slice Data

The cost of obtaining 2D slices is much lower than the cost of directly measuring 3D digital cores. Moreover, rock microfeatures that are more block-like can be obtained through the use of 2D features at the same cost, thereby increasing the flexibility of reconstructing 3D digital cores. Thus, studying 3D digital core-reconstruction methods based on 2D slices is meaningful. This section discusses several challenges that are encountered in this approach.


#### 4.2.1. Importance, Challenges, and Developmental Trends of the Evaluation Function

Two purposes exist for the use of various evaluation functions. One is to implement a stochastic reconstruction method for the statistical constraint of mathematical features; these features can be constrained by an evaluation function to obtain a more accurate reconstruction result. Second, after the digital core is reconstructed, the effectiveness of the reconstruction is tested by using various functions. The ability of the evaluation function to obtain an abstraction of the digital core determines the accuracy of the reconstructed digital core. The evaluation function is also appropriate for determining the number of iterations for the algorithm and the reconstruction speed. The importance of the evaluation function is self-evident. However, our predecessors did not associate an adequate importance with the evaluation function. Therefore, this paper summarizes the evaluation function and provides a reasonable interpretation. Figure15 illustrates the importance of the evaluation function in statistical-based digital core reconstruction.Figure 15
Abstraction of the stochastic reconstruction method based on 2D slice statistical constraint features. The effect of 2D slice to 3D digital core is determined by the evaluation function and algorithm. The evaluation function is one of the important parameters that determine the effect of 3D digital core construction.The variance within a variogram, which is the regionalized variable delta, was proposed by Baniassadi et al. in 1965 [199]. The variogram can be used to test the spatial similarity and internal structural similarity between a reconstructed core and a real core to constrain the reconstructed 3D digital core. The effectiveness of the evaluation is similar to that of the two-point probability function. The function that is used to constrain the porosity of the digital core, which is also known as a single-point probability function, is the most widely used evaluation function. In addition to the single-point probability function, the two-point probability function that was proposed by Smith in 1988 is a classic evaluation function. The probability that any two pixel points are separated by a distance r in a multiphase system is calculated and simultaneously distributed in the same phase to characterize the rock [200]. This function can characterize the spatial distribution of the system. Many scholars considered the use of this function in the establishment of digital cores with different lithologies at different scales [201–204]. In 1992, Lu and Torquato proposed a linear path function that calculates the probability that all voxel points will be in phase with any voxel point at a distance r [205]. This function can reflect the local connectivity (mainly the linear connectivity) of the same phase, which is the same as the two-point probability function, and is a commonly used function in digital core reconstruction [114, 206–210]. The two-point autocorrelation function was proposed in 1995 in a form similar to that of the two-point probability function; both of these functions reflect the in-phase probability of two points at a distance r with only slight differences in the calculation method [211]. Similar to the two-point probability function, the two-point autocorrelation function is also used to reflect the spatial distribution of phases. The N-point probability function was proposed by Roberts in 1997 [212]; this function considers the multiphase conditions of digital cores and extends the single-point and two-point probability functions. The N-point probability function calculates the in-phase probability of each phase that is separated by r, thereby reflecting a higher-order correlation. The N-point probability function contains more details regarding the spatial distribution morphology of the secondary phase. In 1998, a combination of the pore size distribution function and the cumulative pore size distribution function was proposed by Yeong and Torquato [105] to ensure that the pore diameter distribution was reasonable. The cumulative pore size distribution function, which is similar to a linear path function, was based on the proposed pore size distribution function, in which the radius increases and then gradually decreases. In 2000, the local porosity distribution function was proposed by Hilfer [213]; this function calculation method, which can constrain the geometric properties of the reconstructed rock, is more complicated than previously proposed functions. Hajizadeh et al. proposed the connectivity function in 2011 to calculate the pore connectivity of reconstructed and training images. The effectiveness of this function is similar to that of functions that reflect the connectivity of the linear path function. In the same year, Piasecki [214] proposed an entropic descriptor to constrain reconstructed 3D digital cores and improve the speed and accuracy of the 3D digital core reconstruction. In 2014, Ju et al. proposed a fractal system control function based on the general fractal structure of most rocks. In 2015, Karsanina et al. proposed two-point cluster functions to evaluate the reconstruction results [215]; this function first creates a mesh by using the box-covering method and then uses the fractal system control function to describe the fractal. Good results have been obtained when a 3D digital core is reconstructed with this function. Moreover, the fractal system control function, which is different from the previous evaluation function calculation methods, is a new concept. Generally, the evaluation function has more complicated calculation characteristics, and the factors are more carefully considered. The above evaluation functions are organized and shown in Table 4.Table 4
List of evaluation functions of digital cores.
Restricted mathematical properties
Representative evaluation function
Is it common?
Size of the phase
Single point probability function
Yes
Spatial distribution characteristics
Two-point probability function, two-point autocorrelation function,N-point probability function, variogram function
Yes
Phase connectivity
Linear path function, pore size distribution function, cumulative pore size distribution function, connectivity function
Yes
Geometric characteristics
Local porosity function
No
Fractal characteristics
Fractal system control function
NoAlthough the majority of the evaluation functions in Table4 have different formulas, their constrained mathematical properties are very similar. The functions of these evaluation functions overlap with few differences. Hence, evaluation functions can be classified based on the mathematical property that is being constrained: the phase volume size, phase connectivity, spatial distribution characteristics, geometric characteristics, and fractal characteristics. Any similarity between these evaluation functions can increase the evaluation function when reconstructing 3D digital cores. Based on previous research, in which 3D digital cores were reconstructed by the statistic-constrained stochastic reconstruction method, the evaluation function produces results that are very similar to the original 2D slices. This finding indicates that the reconstruction algorithm itself does not have a poor reconstruction effect, but the connectivity of the reconstructed 3D digital core is worse than that of a 3D digital core that is reconstructed through a physical experiment. Finally, existing evaluation functions may not be well suited for the abstraction of shale rocks with extremely complex pore structures and strong heterogeneity. Based on the above problems, some scholars have studied new reconstruction methods; alternatively, stochastic reconstruction methods that are constrained by either the statistics or the morphological characteristics of 2D slices and process-based stochastic reconstruction methods have been used. The above description of the current applicability of the evaluation function is not ideal. In addition, recent studies have been more inclined to employ simulation experiments to test the reconstruction of digital cores, indicating that the ability of the evaluation function to accurately reflect the characteristics of the digital core must be strengthened.The above problem is caused by a poor understanding of the evaluation function in present-day research. First, the existing numerical core evaluation function refers to isotropy; that is, the evaluation function has no directional property, and most evaluation functions are related only to the voxelr. This factor is a serious problem for rocks with complicated pore structures, such as carbonates, igneous rocks, and shales, which display significant heterogeneity, so the effectiveness of the reconstruction is greatly reduced. Second, the correlation among the numerical core evaluation function, pore connectivity, and pore distribution features is presently indirect; this correlation is a voxel-based evaluation method and does not have any petrophysical significance. Therefore, we suggest that the petrophysical characteristics of rocks should be combined (such as in the fractal system control function) and that the directionality or tracking of a certain phase should be considered to achieve a more complex evaluation-function design. This approach should improve the effectiveness of statistically constrained stochastic reconstruction and provide a more accurate assessment of digital cores that are reconstructed by various methods. Equally, the inclusion of more advanced and faster optimization functions, as discussed in Section (4), is also crucial to improving the effectiveness of reconstruction algorithms and developing better mechanisms for reconstructing 3D digital cores. Mathematically and statistically constrained stochastic reconstruction algorithms still have much room for improvement.


#### 4.2.2. Developmental Trends and Challenges of Multiscale 3D Digital Cores

The multiscale problem involves the acquisition of more accurate digital cores with a number of various-resolution images or ultra-high-resolution images for various rocks with different pore types, sizes, and anisotropic characteristics. The pore types of conventional high-porosity and high-permeability sandstones are very simple, thereby eliminating the need to consider multiscale problems. However, the pore origin and pore size vary greatly for unconventional reservoir rocks such as tight sandstones (with a diverse secondary porosity and varying pore radii), carbonates (with large pores such as fissures and voids), and shales (with large differences in the pore sizes between organic and inorganic pores). Therefore, constructing an accurate 3D digital core by using only one type of 2D slice is difficult. For carbonate reservoirs that contain large-scale fractures, coexisting voids, and small-scale primary pores, one must scan large-scale cores with an extremely high resolution to fully reflect the pore structure of the rock; unfortunately, this method is extremely expensive and consumes an excessive amount of memory. When using SEM and other experiments to establish 2D slices, large-scale cracks and holes cannot be observed. Hence, the simulation results of small-scale digital cores are not representative for shale reservoirs because of their heterogeneity. However, the organic pore size of a shale rock is excessively small, so this type of rock is of limited significance for constructing large-scale digital cores. Therefore, the construction of a multiscale digital core is imperative.In 2012, Wang et al. proposed the use of the MCMC method to construct corresponding macroporous digital cores and microporous digital cores and obtained double-peaked carbonate digital cores. In 2013, Yao et al. [216] proposed the reconstruction of 3D digital cores of carbonate reservoirs by using the SA and MCMC methods; macropores were constructed by using the SA method (this method reflects more of the mathematical characteristics, which are more important for macropores, of a 2D slice) and then reconstructed into small pores by using the MCMC method (for small pores, the morphology is more important), after which the pores were finally nested. In 2014, Ma et al. [217] proposed a multiscale digital core-reconstruction method for shale gas rocks. Hebert et al. [218] employed a multiscale method to evaluate the porosity of carbonate rocks. In 2015, Yang et al. [219] proposed a method for combining small-scale digital cores that is equivalent to converting the resulting 3D digital cores into larger-scale 3D digital cores. That same year, Gerke et al. developed a general solution for merging multiscale categorical spatial data into a single dataset by using stochastic reconstructions with rescaled correlation functions [220]. In 2018, Liu et al. [221] developed a multipoint statistical construction method for digital cores that considers microcracks in two steps. Mehmani et al. [222] constructed multiscale digital cores in various diagenetic stages. Sun et al. [223] presented the DRA upscaling method to construct carbonate rock cores by using CT images and slices.Yang et al. [153] postulated that organic pores could be reconstructed by reconstructing inorganic pores with the MPS and MCMC methods; to this end, many organic and inorganic pores were numerically embedded. Figure 16 shows a schematic diagram of a digital core with organic pores, a digital core with inorganic pores, and a nested 3D digital core, while Figure 17 shows a comparison between the pore size distribution curves of the embedded digital cores and the actual experimental pore size distribution curves.Figure 16
Shale 3D digital cores reconstructed using multiscale reconstruction method. Among them, the grey part is the skeleton, and the black part is the pores. (a) Shows an inorganic pore digital core, (b) shows an organic pore digital core, and (c) shows a nested digital core. Figure from [153].
(a)
(b)
(c)Figure 17
Comparing the pore size distribution of the 3D digital core obtained in Figure13 with the experimental pore size distribution, it can be seen that the pore size distribution of the 3D digital core after nesting is in good agreement with the experimental data. Figure reprinted from [153].Research on multiscale issues is still in its infancy, leaving much to be studied and many issues to be discussed. First, a pore size distribution that matches the aperture distribution from a real experiment can be obtained for a numerical core from the superposition nesting method, but the evaluation function and simulation experiments are still inadequate. Questions often arise with regard to how the connectivity of the 3D digital core was reconstructed, how nesting should occur if the simulation results are unsatisfactory, and the reliability of direct nesting, which is equivalent to different holes within the 3D digital core that are not related to its generation. At present, when multiscale digital core reconstruction is performed, the matching and connectivity of different pore sizes are considered relatively low. These concepts are fundamental to matching a digital core with an actual rock, but they are also the most difficult aspect to study, so in-depth research must be conducted. Otherwise, the simulation results of constructed multiscale digital cores will be very inadequate because of the poor connectivity of the pore system, especially for experiments that are related to the pore structure and pore connectivity. We suggest that one must consider whether the obtained results satisfy the spatial relationships between different types of pores and minerals when performing multiscale digital core reconstruction. For example, organic pores in shale rocks are usually spatially similar to organic matter. Cracks that develop in clastic or carbonate rocks usually depend on brittle minerals and have a poor relationship with plastic minerals such as clay. These details can fundamentally improve the stability and reliability of digital cores. Many other details must be considered when reconstructing multiscale digital cores. In summary, one should pay close attention to the matching of different scales when reconstructing multiscale digital cores. These problems are very challenging issues that will be addressed in a future study.


#### 4.2.3. Developmental Trends and Challenges of Multicomponent 3D Digital Cores

The types of minerals in shale reservoir rocks and igneous rocks are diverse, and the mineral, hydrocarbon, and water distributions are complicated; nevertheless, the types and distributions of minerals affect the simulation results, as do the distributions of oil, gas, and water. Thus, one cannot consider only skeleton and pore phases or only skeleton, oil, and water phases for digital cores. In 2011, Jiang et al. [224] added a water-film phase to the reconstructed 3D digital core after the oil-water distribution was constructed; these researchers asserted that the reservoir wettability and natural gas solubility must be considered in the formation water through electrical modelling when evaluating the actual reservoir saturation. In 2016, Nie et al. proposed a multicomponent reconstruction technique that uses the MCMC reconstruction algorithm and nested and obtained a large-scale multicomponent 3D digital core; these researchers also conducted an electrical simulation of the reconstructed 3D digital core and concluded that Archie’s formula is not suitable for shale reservoirs. Zhu and Yu proposed a method for generating cementitious phases within an original reconstructed digital sandstone core. Saad et al. [225] used a physics-based process approach to construct a multicomponent digital core. By simulating a 3D digital core, the cement is considered to influence the moduli of the longitudinal and transverse waves and thus should be considered. Although relatively little attention has been given to multicomponent problems, some scholars have studied the construction of multicomponent 3D digital cores.At present, multicomponent digital core reconstruction has received less attention than multiscale digital core reconstruction because the latter is a common problem that is encountered for many rock types and has yet to be solved. However, many rocks, such as shales, igneous rocks, and metamorphic rocks, typically have multiple components, although not all rocks urgently require this problem to be resolved. The correlations among minerals and between minerals and pores must also be considered to determine whether a reconstructed multicomponent digital core is reasonable and similar to the real core. In other words, the most important component in the construction of multicomponent digital cores is the need to conform to the corresponding microscopic geological features of the rock. Therefore, as with multiscale digital core construction, matching remains the key to the accurate construction of multicomponent digital cores. Similar to the multiscale problem, too few simulation experiments have been performed on multicomponent 3D digital cores; therefore, whether multicomponent 3D digital cores can be directly used in numerical simulation experiments remains unclear. In addition, most existing methods for constructing multicomponent digital cores are directly nested after reconstruction, and the correlation between the geology and minerals is not considered (e.g., organic matter and pyrite are often associated with each other). Therefore, the study of multicomponent reconstruction, which directly uses a reconstruction algorithm, is suggested. In theory, both the MPS method and the SA method can be utilized to reconstruct 3D digital cores that contain more than two phases. In conclusion, multicomponent 3D digital cores are difficult to address and require additional investigation in the future, similarly to multiscale 3D digital cores.


#### 4.2.4. Digital Core Reconstruction Based on Interdisciplinary Research

In addition to the abovementioned issues, this study separately addresses digital core-reconstruction methods based on an interdisciplinary approach to emphasize their importance in the 3D digital core-reconstruction field and discern the potential of interdisciplinary digital core-reconstruction methods.Interdisciplinary research refers to emerging fields in which different disciplines intersect and merge; here, this term mainly refers to the use of other methods to interpret digital core-reconstruction technologies and ideas with the help of other advanced or flourishing methods in combination with the actual need to improve these methods or in combination with classic methods to further improve the reconstruction effect. The most important points when using knowledge and methods of other disciplines for the reconstruction of 3D digital cores are as follows. First, we must pay attention to the progress and trends of other related disciplines, such as image processing, probability theory, modern optimization theory, and machine learning. This approach will be very helpful to further stimulate our thinking and improve the effects of rock reconstruction endeavours. Second, if an interdisciplinary method is required to reconstruct a digital core, one should examine the reconstruction process from an interdisciplinary perspective, which is very crucial. Otherwise, understanding the actual needs of the reconstruction becomes more difficult, leading to inconsistency among the disciplines, as does achieving good results with a new method. Third, one should understand whether the required conditions to use the new method are consistent with the requirements to reconstruct the digital core. If these conditions are not consistent, they must be thoroughly explored and rationally improved with regard to specific problems, which is the key to implementing a new interdisciplinary method.Currently, researching digital core construction has great potential with regard to optimization, particle accumulation, and machine learning. In Section4, we addressed how important the optimization algorithm is for achieving an adequate reconstruction effect; that is, a good optimization algorithm can improve the effectiveness of the reconstruction. An optimization algorithm can also constitute an attempt to improve the accuracy of the 3D digital core reconstruction with interdisciplinary techniques. From the perspective of optimization, the main problems of digital core reconstruction are constraining this process by using several predetermined functions and optimizing the fitness function by changing the values of the data. The digital core that corresponds to the optimal result is the best digital core from a given method. For optimization theory, the digital core-reconstruction problem is only an optimization problem. Therefore, the evaluation function, which determines the reconstruction of the digital core after optimization, becomes very important, which matches the view that we presented in Section 4.2.1.The study of particle accumulation could be very beneficial for digital core-reconstruction methods based on process methods. Many new particle accumulation methods can be directly applied to the reconstruction of digital cores after a simple modification [226–233] because no fundamental difference exists between the accumulation of rock particles and the accumulation of other irregular particles, that is, those deposited by the actions of various forces. Moreover, digital cores are small, so considering influencing factors such as the rhythm and source direction is unnecessary. At present, digital cores that are reconstructed based on the DEM method as proposed by molecular dynamics are more reliable. However, no reports have been published on the use of irregular rock particles for direct DEM-based reconstructions. For instance, Tahmasebi attempted to use only non-DEM irregular particle accumulation methods [234], while other scholars that specialize in particle accumulation have proposed a method for the accumulation of particles with irregular shapes based on DEMs [235]. Thus, the current research on digital core reconstruction is not sufficient.Machine learning is the multidisciplinary and interdisciplinary study of learning, acquiring new knowledge and skills and identifying existing knowledge. Machine learning first attempts to learn the existing rules, after which known rules are employed to predict and judge new samples. Digital core reconstruction could be regarded as such a problem. If the reconstruction law can be learned, then a stable prediction model can be predicted (i.e., a digital core can be reconstructed).Sundararaghavan and Zabaras [236] proposed a method for reconstructing 3D digital cores by using the support vector machine (SVM) method in machine learning. With the continuous development of artificial intelligence and other disciplines, existing algorithms have greatly improved the abilities of pattern recognition and data-reconstruction algorithms compared to the original techniques. Deep learning [237, 238], extreme learning [239], intensive learning [240], and active learning [241] algorithms and generative adversarial networks [242], among others, all boast strong prediction and reconstruction capabilities. The essence of 3D digital core reconstruction is to use known data to obtain the reconstruction law and predict unknown spaces, which coincides with the nature of machine learning algorithms. Currently, the more difficult problem is how to obtain reliable prediction models for high-dimensional and demanding digital core-reconstruction problems with relatively few samples (i.e., when acquiring experimental samples is expensive and time-consuming). The current research on this aspect is still nascent, so evaluating the effect of a small sample size on the actual reconstruction is difficult, constituting another topic to be studied.The above three aspects represent only a very simple discussion of the effects of an interdisciplinary approach on the reconstruction of digital cores. According to the above discussion, a clear interdisciplinary understanding of this problem can significantly improve the effect of digital core reconstruction and expand the application scope of digital core-reconstruction methods (i.e., current digital core-reconstruction methods are relatively simple).Thus, drawing on interdisciplinary content and reasonably utilizing such algorithms is another research direction for 3D digital core reconstruction.


## 4.1. Challenges of Reconstructing 3D Digital Cores Based on Physical Experimental Methods

Existing physical experimental methods to reconstruct 3D digital cores encounter fewer challenges than the creation of 3D digital cores based on 2D slices because the former produce more realistic results. Presently, three key problems must be further studied to establish a 3D digital-core method based on physical experimentation. First, the selected rock sample size and optimal scan resolution is important. For a rock with a small pore size, high-resolution instruments can be used to accurately identify the pore space; conversely, a rock with a larger pore size is recommended for low-resolution pore space identification approaches. Second, binary images require further study. The porosity cannot be determined with a single threshold because of obscure skeleton-porosity boundaries in greyscale images. Additionally, the porosity cannot be determined based on physical property experiments, and the microporosity may be difficult to measure through porosity experiments [196–198]. Therefore, the determination of the segmentation threshold is another problem to be addressed. Finally, the size of the 3D digital core must be modest because of computational limitations. Thus, the determination of an appropriate 3D digital-core size for simulation experiments is also meaningful, and the task of building 3D digital cores based on physical experimentation is far less challenging than the reconstruction of 3D digital cores based on 2D slices.


## 4.2. Trends and Challenges of Reconstructing 3D Digital Cores Based on 2D Slice Data

The cost of obtaining 2D slices is much lower than the cost of directly measuring 3D digital cores. Moreover, rock microfeatures that are more block-like can be obtained through the use of 2D features at the same cost, thereby increasing the flexibility of reconstructing 3D digital cores. Thus, studying 3D digital core-reconstruction methods based on 2D slices is meaningful. This section discusses several challenges that are encountered in this approach.


### 4.2.1. Importance, Challenges, and Developmental Trends of the Evaluation Function

Two purposes exist for the use of various evaluation functions. One is to implement a stochastic reconstruction method for the statistical constraint of mathematical features; these features can be constrained by an evaluation function to obtain a more accurate reconstruction result. Second, after the digital core is reconstructed, the effectiveness of the reconstruction is tested by using various functions. The ability of the evaluation function to obtain an abstraction of the digital core determines the accuracy of the reconstructed digital core. The evaluation function is also appropriate for determining the number of iterations for the algorithm and the reconstruction speed. The importance of the evaluation function is self-evident. However, our predecessors did not associate an adequate importance with the evaluation function. Therefore, this paper summarizes the evaluation function and provides a reasonable interpretation. Figure15 illustrates the importance of the evaluation function in statistical-based digital core reconstruction.Figure 15
Abstraction of the stochastic reconstruction method based on 2D slice statistical constraint features. The effect of 2D slice to 3D digital core is determined by the evaluation function and algorithm. The evaluation function is one of the important parameters that determine the effect of 3D digital core construction.The variance within a variogram, which is the regionalized variable delta, was proposed by Baniassadi et al. in 1965 [199]. The variogram can be used to test the spatial similarity and internal structural similarity between a reconstructed core and a real core to constrain the reconstructed 3D digital core. The effectiveness of the evaluation is similar to that of the two-point probability function. The function that is used to constrain the porosity of the digital core, which is also known as a single-point probability function, is the most widely used evaluation function. In addition to the single-point probability function, the two-point probability function that was proposed by Smith in 1988 is a classic evaluation function. The probability that any two pixel points are separated by a distance r in a multiphase system is calculated and simultaneously distributed in the same phase to characterize the rock [200]. This function can characterize the spatial distribution of the system. Many scholars considered the use of this function in the establishment of digital cores with different lithologies at different scales [201–204]. In 1992, Lu and Torquato proposed a linear path function that calculates the probability that all voxel points will be in phase with any voxel point at a distance r [205]. This function can reflect the local connectivity (mainly the linear connectivity) of the same phase, which is the same as the two-point probability function, and is a commonly used function in digital core reconstruction [114, 206–210]. The two-point autocorrelation function was proposed in 1995 in a form similar to that of the two-point probability function; both of these functions reflect the in-phase probability of two points at a distance r with only slight differences in the calculation method [211]. Similar to the two-point probability function, the two-point autocorrelation function is also used to reflect the spatial distribution of phases. The N-point probability function was proposed by Roberts in 1997 [212]; this function considers the multiphase conditions of digital cores and extends the single-point and two-point probability functions. The N-point probability function calculates the in-phase probability of each phase that is separated by r, thereby reflecting a higher-order correlation. The N-point probability function contains more details regarding the spatial distribution morphology of the secondary phase. In 1998, a combination of the pore size distribution function and the cumulative pore size distribution function was proposed by Yeong and Torquato [105] to ensure that the pore diameter distribution was reasonable. The cumulative pore size distribution function, which is similar to a linear path function, was based on the proposed pore size distribution function, in which the radius increases and then gradually decreases. In 2000, the local porosity distribution function was proposed by Hilfer [213]; this function calculation method, which can constrain the geometric properties of the reconstructed rock, is more complicated than previously proposed functions. Hajizadeh et al. proposed the connectivity function in 2011 to calculate the pore connectivity of reconstructed and training images. The effectiveness of this function is similar to that of functions that reflect the connectivity of the linear path function. In the same year, Piasecki [214] proposed an entropic descriptor to constrain reconstructed 3D digital cores and improve the speed and accuracy of the 3D digital core reconstruction. In 2014, Ju et al. proposed a fractal system control function based on the general fractal structure of most rocks. In 2015, Karsanina et al. proposed two-point cluster functions to evaluate the reconstruction results [215]; this function first creates a mesh by using the box-covering method and then uses the fractal system control function to describe the fractal. Good results have been obtained when a 3D digital core is reconstructed with this function. Moreover, the fractal system control function, which is different from the previous evaluation function calculation methods, is a new concept. Generally, the evaluation function has more complicated calculation characteristics, and the factors are more carefully considered. The above evaluation functions are organized and shown in Table 4.Table 4
List of evaluation functions of digital cores.
Restricted mathematical properties
Representative evaluation function
Is it common?
Size of the phase
Single point probability function
Yes
Spatial distribution characteristics
Two-point probability function, two-point autocorrelation function,N-point probability function, variogram function
Yes
Phase connectivity
Linear path function, pore size distribution function, cumulative pore size distribution function, connectivity function
Yes
Geometric characteristics
Local porosity function
No
Fractal characteristics
Fractal system control function
NoAlthough the majority of the evaluation functions in Table4 have different formulas, their constrained mathematical properties are very similar. The functions of these evaluation functions overlap with few differences. Hence, evaluation functions can be classified based on the mathematical property that is being constrained: the phase volume size, phase connectivity, spatial distribution characteristics, geometric characteristics, and fractal characteristics. Any similarity between these evaluation functions can increase the evaluation function when reconstructing 3D digital cores. Based on previous research, in which 3D digital cores were reconstructed by the statistic-constrained stochastic reconstruction method, the evaluation function produces results that are very similar to the original 2D slices. This finding indicates that the reconstruction algorithm itself does not have a poor reconstruction effect, but the connectivity of the reconstructed 3D digital core is worse than that of a 3D digital core that is reconstructed through a physical experiment. Finally, existing evaluation functions may not be well suited for the abstraction of shale rocks with extremely complex pore structures and strong heterogeneity. Based on the above problems, some scholars have studied new reconstruction methods; alternatively, stochastic reconstruction methods that are constrained by either the statistics or the morphological characteristics of 2D slices and process-based stochastic reconstruction methods have been used. The above description of the current applicability of the evaluation function is not ideal. In addition, recent studies have been more inclined to employ simulation experiments to test the reconstruction of digital cores, indicating that the ability of the evaluation function to accurately reflect the characteristics of the digital core must be strengthened.The above problem is caused by a poor understanding of the evaluation function in present-day research. First, the existing numerical core evaluation function refers to isotropy; that is, the evaluation function has no directional property, and most evaluation functions are related only to the voxelr. This factor is a serious problem for rocks with complicated pore structures, such as carbonates, igneous rocks, and shales, which display significant heterogeneity, so the effectiveness of the reconstruction is greatly reduced. Second, the correlation among the numerical core evaluation function, pore connectivity, and pore distribution features is presently indirect; this correlation is a voxel-based evaluation method and does not have any petrophysical significance. Therefore, we suggest that the petrophysical characteristics of rocks should be combined (such as in the fractal system control function) and that the directionality or tracking of a certain phase should be considered to achieve a more complex evaluation-function design. This approach should improve the effectiveness of statistically constrained stochastic reconstruction and provide a more accurate assessment of digital cores that are reconstructed by various methods. Equally, the inclusion of more advanced and faster optimization functions, as discussed in Section (4), is also crucial to improving the effectiveness of reconstruction algorithms and developing better mechanisms for reconstructing 3D digital cores. Mathematically and statistically constrained stochastic reconstruction algorithms still have much room for improvement.


### 4.2.2. Developmental Trends and Challenges of Multiscale 3D Digital Cores

The multiscale problem involves the acquisition of more accurate digital cores with a number of various-resolution images or ultra-high-resolution images for various rocks with different pore types, sizes, and anisotropic characteristics. The pore types of conventional high-porosity and high-permeability sandstones are very simple, thereby eliminating the need to consider multiscale problems. However, the pore origin and pore size vary greatly for unconventional reservoir rocks such as tight sandstones (with a diverse secondary porosity and varying pore radii), carbonates (with large pores such as fissures and voids), and shales (with large differences in the pore sizes between organic and inorganic pores). Therefore, constructing an accurate 3D digital core by using only one type of 2D slice is difficult. For carbonate reservoirs that contain large-scale fractures, coexisting voids, and small-scale primary pores, one must scan large-scale cores with an extremely high resolution to fully reflect the pore structure of the rock; unfortunately, this method is extremely expensive and consumes an excessive amount of memory. When using SEM and other experiments to establish 2D slices, large-scale cracks and holes cannot be observed. Hence, the simulation results of small-scale digital cores are not representative for shale reservoirs because of their heterogeneity. However, the organic pore size of a shale rock is excessively small, so this type of rock is of limited significance for constructing large-scale digital cores. Therefore, the construction of a multiscale digital core is imperative.In 2012, Wang et al. proposed the use of the MCMC method to construct corresponding macroporous digital cores and microporous digital cores and obtained double-peaked carbonate digital cores. In 2013, Yao et al. [216] proposed the reconstruction of 3D digital cores of carbonate reservoirs by using the SA and MCMC methods; macropores were constructed by using the SA method (this method reflects more of the mathematical characteristics, which are more important for macropores, of a 2D slice) and then reconstructed into small pores by using the MCMC method (for small pores, the morphology is more important), after which the pores were finally nested. In 2014, Ma et al. [217] proposed a multiscale digital core-reconstruction method for shale gas rocks. Hebert et al. [218] employed a multiscale method to evaluate the porosity of carbonate rocks. In 2015, Yang et al. [219] proposed a method for combining small-scale digital cores that is equivalent to converting the resulting 3D digital cores into larger-scale 3D digital cores. That same year, Gerke et al. developed a general solution for merging multiscale categorical spatial data into a single dataset by using stochastic reconstructions with rescaled correlation functions [220]. In 2018, Liu et al. [221] developed a multipoint statistical construction method for digital cores that considers microcracks in two steps. Mehmani et al. [222] constructed multiscale digital cores in various diagenetic stages. Sun et al. [223] presented the DRA upscaling method to construct carbonate rock cores by using CT images and slices.Yang et al. [153] postulated that organic pores could be reconstructed by reconstructing inorganic pores with the MPS and MCMC methods; to this end, many organic and inorganic pores were numerically embedded. Figure 16 shows a schematic diagram of a digital core with organic pores, a digital core with inorganic pores, and a nested 3D digital core, while Figure 17 shows a comparison between the pore size distribution curves of the embedded digital cores and the actual experimental pore size distribution curves.Figure 16
Shale 3D digital cores reconstructed using multiscale reconstruction method. Among them, the grey part is the skeleton, and the black part is the pores. (a) Shows an inorganic pore digital core, (b) shows an organic pore digital core, and (c) shows a nested digital core. Figure from [153].
(a)
(b)
(c)Figure 17
Comparing the pore size distribution of the 3D digital core obtained in Figure13 with the experimental pore size distribution, it can be seen that the pore size distribution of the 3D digital core after nesting is in good agreement with the experimental data. Figure reprinted from [153].Research on multiscale issues is still in its infancy, leaving much to be studied and many issues to be discussed. First, a pore size distribution that matches the aperture distribution from a real experiment can be obtained for a numerical core from the superposition nesting method, but the evaluation function and simulation experiments are still inadequate. Questions often arise with regard to how the connectivity of the 3D digital core was reconstructed, how nesting should occur if the simulation results are unsatisfactory, and the reliability of direct nesting, which is equivalent to different holes within the 3D digital core that are not related to its generation. At present, when multiscale digital core reconstruction is performed, the matching and connectivity of different pore sizes are considered relatively low. These concepts are fundamental to matching a digital core with an actual rock, but they are also the most difficult aspect to study, so in-depth research must be conducted. Otherwise, the simulation results of constructed multiscale digital cores will be very inadequate because of the poor connectivity of the pore system, especially for experiments that are related to the pore structure and pore connectivity. We suggest that one must consider whether the obtained results satisfy the spatial relationships between different types of pores and minerals when performing multiscale digital core reconstruction. For example, organic pores in shale rocks are usually spatially similar to organic matter. Cracks that develop in clastic or carbonate rocks usually depend on brittle minerals and have a poor relationship with plastic minerals such as clay. These details can fundamentally improve the stability and reliability of digital cores. Many other details must be considered when reconstructing multiscale digital cores. In summary, one should pay close attention to the matching of different scales when reconstructing multiscale digital cores. These problems are very challenging issues that will be addressed in a future study.


### 4.2.3. Developmental Trends and Challenges of Multicomponent 3D Digital Cores

The types of minerals in shale reservoir rocks and igneous rocks are diverse, and the mineral, hydrocarbon, and water distributions are complicated; nevertheless, the types and distributions of minerals affect the simulation results, as do the distributions of oil, gas, and water. Thus, one cannot consider only skeleton and pore phases or only skeleton, oil, and water phases for digital cores. In 2011, Jiang et al. [224] added a water-film phase to the reconstructed 3D digital core after the oil-water distribution was constructed; these researchers asserted that the reservoir wettability and natural gas solubility must be considered in the formation water through electrical modelling when evaluating the actual reservoir saturation. In 2016, Nie et al. proposed a multicomponent reconstruction technique that uses the MCMC reconstruction algorithm and nested and obtained a large-scale multicomponent 3D digital core; these researchers also conducted an electrical simulation of the reconstructed 3D digital core and concluded that Archie’s formula is not suitable for shale reservoirs. Zhu and Yu proposed a method for generating cementitious phases within an original reconstructed digital sandstone core. Saad et al. [225] used a physics-based process approach to construct a multicomponent digital core. By simulating a 3D digital core, the cement is considered to influence the moduli of the longitudinal and transverse waves and thus should be considered. Although relatively little attention has been given to multicomponent problems, some scholars have studied the construction of multicomponent 3D digital cores.At present, multicomponent digital core reconstruction has received less attention than multiscale digital core reconstruction because the latter is a common problem that is encountered for many rock types and has yet to be solved. However, many rocks, such as shales, igneous rocks, and metamorphic rocks, typically have multiple components, although not all rocks urgently require this problem to be resolved. The correlations among minerals and between minerals and pores must also be considered to determine whether a reconstructed multicomponent digital core is reasonable and similar to the real core. In other words, the most important component in the construction of multicomponent digital cores is the need to conform to the corresponding microscopic geological features of the rock. Therefore, as with multiscale digital core construction, matching remains the key to the accurate construction of multicomponent digital cores. Similar to the multiscale problem, too few simulation experiments have been performed on multicomponent 3D digital cores; therefore, whether multicomponent 3D digital cores can be directly used in numerical simulation experiments remains unclear. In addition, most existing methods for constructing multicomponent digital cores are directly nested after reconstruction, and the correlation between the geology and minerals is not considered (e.g., organic matter and pyrite are often associated with each other). Therefore, the study of multicomponent reconstruction, which directly uses a reconstruction algorithm, is suggested. In theory, both the MPS method and the SA method can be utilized to reconstruct 3D digital cores that contain more than two phases. In conclusion, multicomponent 3D digital cores are difficult to address and require additional investigation in the future, similarly to multiscale 3D digital cores.


### 4.2.4. Digital Core Reconstruction Based on Interdisciplinary Research

In addition to the abovementioned issues, this study separately addresses digital core-reconstruction methods based on an interdisciplinary approach to emphasize their importance in the 3D digital core-reconstruction field and discern the potential of interdisciplinary digital core-reconstruction methods.Interdisciplinary research refers to emerging fields in which different disciplines intersect and merge; here, this term mainly refers to the use of other methods to interpret digital core-reconstruction technologies and ideas with the help of other advanced or flourishing methods in combination with the actual need to improve these methods or in combination with classic methods to further improve the reconstruction effect. The most important points when using knowledge and methods of other disciplines for the reconstruction of 3D digital cores are as follows. First, we must pay attention to the progress and trends of other related disciplines, such as image processing, probability theory, modern optimization theory, and machine learning. This approach will be very helpful to further stimulate our thinking and improve the effects of rock reconstruction endeavours. Second, if an interdisciplinary method is required to reconstruct a digital core, one should examine the reconstruction process from an interdisciplinary perspective, which is very crucial. Otherwise, understanding the actual needs of the reconstruction becomes more difficult, leading to inconsistency among the disciplines, as does achieving good results with a new method. Third, one should understand whether the required conditions to use the new method are consistent with the requirements to reconstruct the digital core. If these conditions are not consistent, they must be thoroughly explored and rationally improved with regard to specific problems, which is the key to implementing a new interdisciplinary method.Currently, researching digital core construction has great potential with regard to optimization, particle accumulation, and machine learning. In Section4, we addressed how important the optimization algorithm is for achieving an adequate reconstruction effect; that is, a good optimization algorithm can improve the effectiveness of the reconstruction. An optimization algorithm can also constitute an attempt to improve the accuracy of the 3D digital core reconstruction with interdisciplinary techniques. From the perspective of optimization, the main problems of digital core reconstruction are constraining this process by using several predetermined functions and optimizing the fitness function by changing the values of the data. The digital core that corresponds to the optimal result is the best digital core from a given method. For optimization theory, the digital core-reconstruction problem is only an optimization problem. Therefore, the evaluation function, which determines the reconstruction of the digital core after optimization, becomes very important, which matches the view that we presented in Section 4.2.1.The study of particle accumulation could be very beneficial for digital core-reconstruction methods based on process methods. Many new particle accumulation methods can be directly applied to the reconstruction of digital cores after a simple modification [226–233] because no fundamental difference exists between the accumulation of rock particles and the accumulation of other irregular particles, that is, those deposited by the actions of various forces. Moreover, digital cores are small, so considering influencing factors such as the rhythm and source direction is unnecessary. At present, digital cores that are reconstructed based on the DEM method as proposed by molecular dynamics are more reliable. However, no reports have been published on the use of irregular rock particles for direct DEM-based reconstructions. For instance, Tahmasebi attempted to use only non-DEM irregular particle accumulation methods [234], while other scholars that specialize in particle accumulation have proposed a method for the accumulation of particles with irregular shapes based on DEMs [235]. Thus, the current research on digital core reconstruction is not sufficient.Machine learning is the multidisciplinary and interdisciplinary study of learning, acquiring new knowledge and skills and identifying existing knowledge. Machine learning first attempts to learn the existing rules, after which known rules are employed to predict and judge new samples. Digital core reconstruction could be regarded as such a problem. If the reconstruction law can be learned, then a stable prediction model can be predicted (i.e., a digital core can be reconstructed).Sundararaghavan and Zabaras [236] proposed a method for reconstructing 3D digital cores by using the support vector machine (SVM) method in machine learning. With the continuous development of artificial intelligence and other disciplines, existing algorithms have greatly improved the abilities of pattern recognition and data-reconstruction algorithms compared to the original techniques. Deep learning [237, 238], extreme learning [239], intensive learning [240], and active learning [241] algorithms and generative adversarial networks [242], among others, all boast strong prediction and reconstruction capabilities. The essence of 3D digital core reconstruction is to use known data to obtain the reconstruction law and predict unknown spaces, which coincides with the nature of machine learning algorithms. Currently, the more difficult problem is how to obtain reliable prediction models for high-dimensional and demanding digital core-reconstruction problems with relatively few samples (i.e., when acquiring experimental samples is expensive and time-consuming). The current research on this aspect is still nascent, so evaluating the effect of a small sample size on the actual reconstruction is difficult, constituting another topic to be studied.The above three aspects represent only a very simple discussion of the effects of an interdisciplinary approach on the reconstruction of digital cores. According to the above discussion, a clear interdisciplinary understanding of this problem can significantly improve the effect of digital core reconstruction and expand the application scope of digital core-reconstruction methods (i.e., current digital core-reconstruction methods are relatively simple).Thus, drawing on interdisciplinary content and reasonably utilizing such algorithms is another research direction for 3D digital core reconstruction.


## 4.2.1. Importance, Challenges, and Developmental Trends of the Evaluation Function

Two purposes exist for the use of various evaluation functions. One is to implement a stochastic reconstruction method for the statistical constraint of mathematical features; these features can be constrained by an evaluation function to obtain a more accurate reconstruction result. Second, after the digital core is reconstructed, the effectiveness of the reconstruction is tested by using various functions. The ability of the evaluation function to obtain an abstraction of the digital core determines the accuracy of the reconstructed digital core. The evaluation function is also appropriate for determining the number of iterations for the algorithm and the reconstruction speed. The importance of the evaluation function is self-evident. However, our predecessors did not associate an adequate importance with the evaluation function. Therefore, this paper summarizes the evaluation function and provides a reasonable interpretation. Figure15 illustrates the importance of the evaluation function in statistical-based digital core reconstruction.Figure 15
Abstraction of the stochastic reconstruction method based on 2D slice statistical constraint features. The effect of 2D slice to 3D digital core is determined by the evaluation function and algorithm. The evaluation function is one of the important parameters that determine the effect of 3D digital core construction.The variance within a variogram, which is the regionalized variable delta, was proposed by Baniassadi et al. in 1965 [199]. The variogram can be used to test the spatial similarity and internal structural similarity between a reconstructed core and a real core to constrain the reconstructed 3D digital core. The effectiveness of the evaluation is similar to that of the two-point probability function. The function that is used to constrain the porosity of the digital core, which is also known as a single-point probability function, is the most widely used evaluation function. In addition to the single-point probability function, the two-point probability function that was proposed by Smith in 1988 is a classic evaluation function. The probability that any two pixel points are separated by a distance r in a multiphase system is calculated and simultaneously distributed in the same phase to characterize the rock [200]. This function can characterize the spatial distribution of the system. Many scholars considered the use of this function in the establishment of digital cores with different lithologies at different scales [201–204]. In 1992, Lu and Torquato proposed a linear path function that calculates the probability that all voxel points will be in phase with any voxel point at a distance r [205]. This function can reflect the local connectivity (mainly the linear connectivity) of the same phase, which is the same as the two-point probability function, and is a commonly used function in digital core reconstruction [114, 206–210]. The two-point autocorrelation function was proposed in 1995 in a form similar to that of the two-point probability function; both of these functions reflect the in-phase probability of two points at a distance r with only slight differences in the calculation method [211]. Similar to the two-point probability function, the two-point autocorrelation function is also used to reflect the spatial distribution of phases. The N-point probability function was proposed by Roberts in 1997 [212]; this function considers the multiphase conditions of digital cores and extends the single-point and two-point probability functions. The N-point probability function calculates the in-phase probability of each phase that is separated by r, thereby reflecting a higher-order correlation. The N-point probability function contains more details regarding the spatial distribution morphology of the secondary phase. In 1998, a combination of the pore size distribution function and the cumulative pore size distribution function was proposed by Yeong and Torquato [105] to ensure that the pore diameter distribution was reasonable. The cumulative pore size distribution function, which is similar to a linear path function, was based on the proposed pore size distribution function, in which the radius increases and then gradually decreases. In 2000, the local porosity distribution function was proposed by Hilfer [213]; this function calculation method, which can constrain the geometric properties of the reconstructed rock, is more complicated than previously proposed functions. Hajizadeh et al. proposed the connectivity function in 2011 to calculate the pore connectivity of reconstructed and training images. The effectiveness of this function is similar to that of functions that reflect the connectivity of the linear path function. In the same year, Piasecki [214] proposed an entropic descriptor to constrain reconstructed 3D digital cores and improve the speed and accuracy of the 3D digital core reconstruction. In 2014, Ju et al. proposed a fractal system control function based on the general fractal structure of most rocks. In 2015, Karsanina et al. proposed two-point cluster functions to evaluate the reconstruction results [215]; this function first creates a mesh by using the box-covering method and then uses the fractal system control function to describe the fractal. Good results have been obtained when a 3D digital core is reconstructed with this function. Moreover, the fractal system control function, which is different from the previous evaluation function calculation methods, is a new concept. Generally, the evaluation function has more complicated calculation characteristics, and the factors are more carefully considered. The above evaluation functions are organized and shown in Table 4.Table 4
List of evaluation functions of digital cores.
Restricted mathematical properties
Representative evaluation function
Is it common?
Size of the phase
Single point probability function
Yes
Spatial distribution characteristics
Two-point probability function, two-point autocorrelation function,N-point probability function, variogram function
Yes
Phase connectivity
Linear path function, pore size distribution function, cumulative pore size distribution function, connectivity function
Yes
Geometric characteristics
Local porosity function
No
Fractal characteristics
Fractal system control function
NoAlthough the majority of the evaluation functions in Table4 have different formulas, their constrained mathematical properties are very similar. The functions of these evaluation functions overlap with few differences. Hence, evaluation functions can be classified based on the mathematical property that is being constrained: the phase volume size, phase connectivity, spatial distribution characteristics, geometric characteristics, and fractal characteristics. Any similarity between these evaluation functions can increase the evaluation function when reconstructing 3D digital cores. Based on previous research, in which 3D digital cores were reconstructed by the statistic-constrained stochastic reconstruction method, the evaluation function produces results that are very similar to the original 2D slices. This finding indicates that the reconstruction algorithm itself does not have a poor reconstruction effect, but the connectivity of the reconstructed 3D digital core is worse than that of a 3D digital core that is reconstructed through a physical experiment. Finally, existing evaluation functions may not be well suited for the abstraction of shale rocks with extremely complex pore structures and strong heterogeneity. Based on the above problems, some scholars have studied new reconstruction methods; alternatively, stochastic reconstruction methods that are constrained by either the statistics or the morphological characteristics of 2D slices and process-based stochastic reconstruction methods have been used. The above description of the current applicability of the evaluation function is not ideal. In addition, recent studies have been more inclined to employ simulation experiments to test the reconstruction of digital cores, indicating that the ability of the evaluation function to accurately reflect the characteristics of the digital core must be strengthened.The above problem is caused by a poor understanding of the evaluation function in present-day research. First, the existing numerical core evaluation function refers to isotropy; that is, the evaluation function has no directional property, and most evaluation functions are related only to the voxelr. This factor is a serious problem for rocks with complicated pore structures, such as carbonates, igneous rocks, and shales, which display significant heterogeneity, so the effectiveness of the reconstruction is greatly reduced. Second, the correlation among the numerical core evaluation function, pore connectivity, and pore distribution features is presently indirect; this correlation is a voxel-based evaluation method and does not have any petrophysical significance. Therefore, we suggest that the petrophysical characteristics of rocks should be combined (such as in the fractal system control function) and that the directionality or tracking of a certain phase should be considered to achieve a more complex evaluation-function design. This approach should improve the effectiveness of statistically constrained stochastic reconstruction and provide a more accurate assessment of digital cores that are reconstructed by various methods. Equally, the inclusion of more advanced and faster optimization functions, as discussed in Section (4), is also crucial to improving the effectiveness of reconstruction algorithms and developing better mechanisms for reconstructing 3D digital cores. Mathematically and statistically constrained stochastic reconstruction algorithms still have much room for improvement.


## 4.2.2. Developmental Trends and Challenges of Multiscale 3D Digital Cores

The multiscale problem involves the acquisition of more accurate digital cores with a number of various-resolution images or ultra-high-resolution images for various rocks with different pore types, sizes, and anisotropic characteristics. The pore types of conventional high-porosity and high-permeability sandstones are very simple, thereby eliminating the need to consider multiscale problems. However, the pore origin and pore size vary greatly for unconventional reservoir rocks such as tight sandstones (with a diverse secondary porosity and varying pore radii), carbonates (with large pores such as fissures and voids), and shales (with large differences in the pore sizes between organic and inorganic pores). Therefore, constructing an accurate 3D digital core by using only one type of 2D slice is difficult. For carbonate reservoirs that contain large-scale fractures, coexisting voids, and small-scale primary pores, one must scan large-scale cores with an extremely high resolution to fully reflect the pore structure of the rock; unfortunately, this method is extremely expensive and consumes an excessive amount of memory. When using SEM and other experiments to establish 2D slices, large-scale cracks and holes cannot be observed. Hence, the simulation results of small-scale digital cores are not representative for shale reservoirs because of their heterogeneity. However, the organic pore size of a shale rock is excessively small, so this type of rock is of limited significance for constructing large-scale digital cores. Therefore, the construction of a multiscale digital core is imperative.In 2012, Wang et al. proposed the use of the MCMC method to construct corresponding macroporous digital cores and microporous digital cores and obtained double-peaked carbonate digital cores. In 2013, Yao et al. [216] proposed the reconstruction of 3D digital cores of carbonate reservoirs by using the SA and MCMC methods; macropores were constructed by using the SA method (this method reflects more of the mathematical characteristics, which are more important for macropores, of a 2D slice) and then reconstructed into small pores by using the MCMC method (for small pores, the morphology is more important), after which the pores were finally nested. In 2014, Ma et al. [217] proposed a multiscale digital core-reconstruction method for shale gas rocks. Hebert et al. [218] employed a multiscale method to evaluate the porosity of carbonate rocks. In 2015, Yang et al. [219] proposed a method for combining small-scale digital cores that is equivalent to converting the resulting 3D digital cores into larger-scale 3D digital cores. That same year, Gerke et al. developed a general solution for merging multiscale categorical spatial data into a single dataset by using stochastic reconstructions with rescaled correlation functions [220]. In 2018, Liu et al. [221] developed a multipoint statistical construction method for digital cores that considers microcracks in two steps. Mehmani et al. [222] constructed multiscale digital cores in various diagenetic stages. Sun et al. [223] presented the DRA upscaling method to construct carbonate rock cores by using CT images and slices.Yang et al. [153] postulated that organic pores could be reconstructed by reconstructing inorganic pores with the MPS and MCMC methods; to this end, many organic and inorganic pores were numerically embedded. Figure 16 shows a schematic diagram of a digital core with organic pores, a digital core with inorganic pores, and a nested 3D digital core, while Figure 17 shows a comparison between the pore size distribution curves of the embedded digital cores and the actual experimental pore size distribution curves.Figure 16
Shale 3D digital cores reconstructed using multiscale reconstruction method. Among them, the grey part is the skeleton, and the black part is the pores. (a) Shows an inorganic pore digital core, (b) shows an organic pore digital core, and (c) shows a nested digital core. Figure from [153].
(a)
(b)
(c)Figure 17
Comparing the pore size distribution of the 3D digital core obtained in Figure13 with the experimental pore size distribution, it can be seen that the pore size distribution of the 3D digital core after nesting is in good agreement with the experimental data. Figure reprinted from [153].Research on multiscale issues is still in its infancy, leaving much to be studied and many issues to be discussed. First, a pore size distribution that matches the aperture distribution from a real experiment can be obtained for a numerical core from the superposition nesting method, but the evaluation function and simulation experiments are still inadequate. Questions often arise with regard to how the connectivity of the 3D digital core was reconstructed, how nesting should occur if the simulation results are unsatisfactory, and the reliability of direct nesting, which is equivalent to different holes within the 3D digital core that are not related to its generation. At present, when multiscale digital core reconstruction is performed, the matching and connectivity of different pore sizes are considered relatively low. These concepts are fundamental to matching a digital core with an actual rock, but they are also the most difficult aspect to study, so in-depth research must be conducted. Otherwise, the simulation results of constructed multiscale digital cores will be very inadequate because of the poor connectivity of the pore system, especially for experiments that are related to the pore structure and pore connectivity. We suggest that one must consider whether the obtained results satisfy the spatial relationships between different types of pores and minerals when performing multiscale digital core reconstruction. For example, organic pores in shale rocks are usually spatially similar to organic matter. Cracks that develop in clastic or carbonate rocks usually depend on brittle minerals and have a poor relationship with plastic minerals such as clay. These details can fundamentally improve the stability and reliability of digital cores. Many other details must be considered when reconstructing multiscale digital cores. In summary, one should pay close attention to the matching of different scales when reconstructing multiscale digital cores. These problems are very challenging issues that will be addressed in a future study.


## 4.2.3. Developmental Trends and Challenges of Multicomponent 3D Digital Cores

The types of minerals in shale reservoir rocks and igneous rocks are diverse, and the mineral, hydrocarbon, and water distributions are complicated; nevertheless, the types and distributions of minerals affect the simulation results, as do the distributions of oil, gas, and water. Thus, one cannot consider only skeleton and pore phases or only skeleton, oil, and water phases for digital cores. In 2011, Jiang et al. [224] added a water-film phase to the reconstructed 3D digital core after the oil-water distribution was constructed; these researchers asserted that the reservoir wettability and natural gas solubility must be considered in the formation water through electrical modelling when evaluating the actual reservoir saturation. In 2016, Nie et al. proposed a multicomponent reconstruction technique that uses the MCMC reconstruction algorithm and nested and obtained a large-scale multicomponent 3D digital core; these researchers also conducted an electrical simulation of the reconstructed 3D digital core and concluded that Archie’s formula is not suitable for shale reservoirs. Zhu and Yu proposed a method for generating cementitious phases within an original reconstructed digital sandstone core. Saad et al. [225] used a physics-based process approach to construct a multicomponent digital core. By simulating a 3D digital core, the cement is considered to influence the moduli of the longitudinal and transverse waves and thus should be considered. Although relatively little attention has been given to multicomponent problems, some scholars have studied the construction of multicomponent 3D digital cores.At present, multicomponent digital core reconstruction has received less attention than multiscale digital core reconstruction because the latter is a common problem that is encountered for many rock types and has yet to be solved. However, many rocks, such as shales, igneous rocks, and metamorphic rocks, typically have multiple components, although not all rocks urgently require this problem to be resolved. The correlations among minerals and between minerals and pores must also be considered to determine whether a reconstructed multicomponent digital core is reasonable and similar to the real core. In other words, the most important component in the construction of multicomponent digital cores is the need to conform to the corresponding microscopic geological features of the rock. Therefore, as with multiscale digital core construction, matching remains the key to the accurate construction of multicomponent digital cores. Similar to the multiscale problem, too few simulation experiments have been performed on multicomponent 3D digital cores; therefore, whether multicomponent 3D digital cores can be directly used in numerical simulation experiments remains unclear. In addition, most existing methods for constructing multicomponent digital cores are directly nested after reconstruction, and the correlation between the geology and minerals is not considered (e.g., organic matter and pyrite are often associated with each other). Therefore, the study of multicomponent reconstruction, which directly uses a reconstruction algorithm, is suggested. In theory, both the MPS method and the SA method can be utilized to reconstruct 3D digital cores that contain more than two phases. In conclusion, multicomponent 3D digital cores are difficult to address and require additional investigation in the future, similarly to multiscale 3D digital cores.


## 4.2.4. Digital Core Reconstruction Based on Interdisciplinary Research

In addition to the abovementioned issues, this study separately addresses digital core-reconstruction methods based on an interdisciplinary approach to emphasize their importance in the 3D digital core-reconstruction field and discern the potential of interdisciplinary digital core-reconstruction methods.Interdisciplinary research refers to emerging fields in which different disciplines intersect and merge; here, this term mainly refers to the use of other methods to interpret digital core-reconstruction technologies and ideas with the help of other advanced or flourishing methods in combination with the actual need to improve these methods or in combination with classic methods to further improve the reconstruction effect. The most important points when using knowledge and methods of other disciplines for the reconstruction of 3D digital cores are as follows. First, we must pay attention to the progress and trends of other related disciplines, such as image processing, probability theory, modern optimization theory, and machine learning. This approach will be very helpful to further stimulate our thinking and improve the effects of rock reconstruction endeavours. Second, if an interdisciplinary method is required to reconstruct a digital core, one should examine the reconstruction process from an interdisciplinary perspective, which is very crucial. Otherwise, understanding the actual needs of the reconstruction becomes more difficult, leading to inconsistency among the disciplines, as does achieving good results with a new method. Third, one should understand whether the required conditions to use the new method are consistent with the requirements to reconstruct the digital core. If these conditions are not consistent, they must be thoroughly explored and rationally improved with regard to specific problems, which is the key to implementing a new interdisciplinary method.Currently, researching digital core construction has great potential with regard to optimization, particle accumulation, and machine learning. In Section4, we addressed how important the optimization algorithm is for achieving an adequate reconstruction effect; that is, a good optimization algorithm can improve the effectiveness of the reconstruction. An optimization algorithm can also constitute an attempt to improve the accuracy of the 3D digital core reconstruction with interdisciplinary techniques. From the perspective of optimization, the main problems of digital core reconstruction are constraining this process by using several predetermined functions and optimizing the fitness function by changing the values of the data. The digital core that corresponds to the optimal result is the best digital core from a given method. For optimization theory, the digital core-reconstruction problem is only an optimization problem. Therefore, the evaluation function, which determines the reconstruction of the digital core after optimization, becomes very important, which matches the view that we presented in Section 4.2.1.The study of particle accumulation could be very beneficial for digital core-reconstruction methods based on process methods. Many new particle accumulation methods can be directly applied to the reconstruction of digital cores after a simple modification [226–233] because no fundamental difference exists between the accumulation of rock particles and the accumulation of other irregular particles, that is, those deposited by the actions of various forces. Moreover, digital cores are small, so considering influencing factors such as the rhythm and source direction is unnecessary. At present, digital cores that are reconstructed based on the DEM method as proposed by molecular dynamics are more reliable. However, no reports have been published on the use of irregular rock particles for direct DEM-based reconstructions. For instance, Tahmasebi attempted to use only non-DEM irregular particle accumulation methods [234], while other scholars that specialize in particle accumulation have proposed a method for the accumulation of particles with irregular shapes based on DEMs [235]. Thus, the current research on digital core reconstruction is not sufficient.Machine learning is the multidisciplinary and interdisciplinary study of learning, acquiring new knowledge and skills and identifying existing knowledge. Machine learning first attempts to learn the existing rules, after which known rules are employed to predict and judge new samples. Digital core reconstruction could be regarded as such a problem. If the reconstruction law can be learned, then a stable prediction model can be predicted (i.e., a digital core can be reconstructed).Sundararaghavan and Zabaras [236] proposed a method for reconstructing 3D digital cores by using the support vector machine (SVM) method in machine learning. With the continuous development of artificial intelligence and other disciplines, existing algorithms have greatly improved the abilities of pattern recognition and data-reconstruction algorithms compared to the original techniques. Deep learning [237, 238], extreme learning [239], intensive learning [240], and active learning [241] algorithms and generative adversarial networks [242], among others, all boast strong prediction and reconstruction capabilities. The essence of 3D digital core reconstruction is to use known data to obtain the reconstruction law and predict unknown spaces, which coincides with the nature of machine learning algorithms. Currently, the more difficult problem is how to obtain reliable prediction models for high-dimensional and demanding digital core-reconstruction problems with relatively few samples (i.e., when acquiring experimental samples is expensive and time-consuming). The current research on this aspect is still nascent, so evaluating the effect of a small sample size on the actual reconstruction is difficult, constituting another topic to be studied.The above three aspects represent only a very simple discussion of the effects of an interdisciplinary approach on the reconstruction of digital cores. According to the above discussion, a clear interdisciplinary understanding of this problem can significantly improve the effect of digital core reconstruction and expand the application scope of digital core-reconstruction methods (i.e., current digital core-reconstruction methods are relatively simple).Thus, drawing on interdisciplinary content and reasonably utilizing such algorithms is another research direction for 3D digital core reconstruction.


## 5. Conclusions

This study provided an overview of 3D digital core-reconstruction methods and investigated the most common techniques that are employed for the reconstruction of 3D digital cores.(1)
Reconstruction methods for digital cores can be divided into two approaches for establishing 3D digital cores: those based on physical experiments and those based on 2D slices. The methods of establishing 3D digital cores based on physical experiments include confocal laser scanning, the imaging method based on a serial section rock overlay, and X-ray CT scanning; among these methods, the highest resolution can reach the nanometre level. The methods for reconstructing 3D digital cores based on 2D slices include statistically constrained stochastic reconstruction based on mathematical features, the statistically constrained stochastic reconstruction of sliced morphological features, process-based stochastic reconstruction, and hybrid stochastic reconstruction
(2)
Common statistical feature-constrained stochastic reconstruction methods include the Gaussian field method, random search algorithm based on multiple operators, SA algorithm, and sequential-indicator simulation technique. Digital cores that are reconstructed with the SA method boast good connectivity and are easy to improve; therefore, this method has the most potential and is very suitable for forming combinations with other types of methods
(3)
The most common stochastic reconstruction methods based on the statistical constraints of 2D slice morphological characteristics are the MCMC reconstruction method and the MPS reconstruction method. These methods can better reproduce the morphology of each 2D slice phase and represent the best reconstruction methods that are currently available. However, these techniques do not consider the mathematical statistics of 2D slices, and some of the statistical information may not effectively reflect the characteristics of the 2D slices. In addition, compared to the statistically constrained stochastic reconstruction method of mathematical features, this method is more difficult to improve
(4)
The connectivity of 3D digital cores from the process-based stochastic reconstruction approach is better than that from other methods. One disadvantage of this method is that the diagenesis process is considered too simple. In the future, the diagenetic processes of different lithologies should be integrated to study the 3D digital core reconstruction of unconventional reservoir rocks
(5)
Hybrid stochastic reconstruction methods combine the advantages of various methods to complement each other and can consider more reconstruction factors than can any single method; thus, its reconstruction speed and reconstruction effect are superior. Hybrid stochastic reconstruction methods have some significance, and improvements to their accuracy depend on both the improvement of a single reconstruction method and the combination of different reconstruction methods
(6)
In terms of the challenges that are encountered during the reconstruction of 3D digital cores, the current challenge of establishing 3D digital core methods based on physical experiments is less pressing than that of establishing 3D digital core methods based on 2D slices. The challenges that are associated with the former are determined by the sample size, the optimal scan rate, the binarization of greyscale images, and the size of the 3D digital core that is used for the simulation. The challenges that are associated with the latter are related to four major aspects: finding a better evaluation function, establishing multiscale digital cores, establishing multicomponent digital cores, and using an interdisciplinary approach to improve the accuracy of 3D digital core reconstruction. Presently, the methods and theories for establishing 3D digital cores based on 2D slices are immature, and many challenges remain. Moreover, the problems that are considered in this research are still relatively simple. Relevant research on 3D digital core-reconstruction methods will have to utilize broader thinking and creative research



---

*Source: 7814180-2019-05-20.xml*