Microscopic Stress Sensitivity Analysis with In Situ SEM Study and Digital Core Deformation Simulation

Rock stress sensitivity is typically investigated macroscopically. In contrast, a new method combining in situ Scanning Electronic Microscope (SEM) study and digital core deformation simulation is developed in this paper, providing an effective way to investigate the relationship between microstructural deformation and decreasing permeability. The simulation method might replace in situ SEM study under certain scenarios. First, the in situ SEM study was implemented, and the microstructure deformations of rock samples during uniaxial loading were observed and recorded. The SEM images at different stress states were analyzed by digital image correlation (DIC) technique to investigate the principles of these deformations. A deformation simulation method was correspondingly proposed. The simulation effectiveness was demonstrated by comparing the simulation and the in situ SEM study results. To validate the simulation method for the three-dimensional (3D) digital core, porosity-permeability integrated measurements under triaxial stresses were conducted to obtain macroscale data under different stress states for a tight sandstone sample. A 3D digital core was reconstructed by micro-CT imaging with the same rock sample. Under the constraints of the measured porosity changes, the 3D digital core deformation was simulated. A series of simulated cores at different stress states were used for pore network model extraction, and the corresponding permeability was calculated. A comparison of the permeability changes of the simulation and porosity-permeability integrated measurements indicated consistently that the simulation method can characterize the 3D digital core stress sensitivity. In addition, the in situ SEM study results revealed that the throats deformed more severely than the pores by generating the pore and throat diameter frequency distributions at different stress states. Therefore, we concluded that throat deformation is more critical than pore deformation for permeability reduction.


Introduction
For the production of unconventional oil and gas in tight sandstone and shale formations, people are mostly concerned with the stress sensitivity effect on reservoir rock permeability [1][2][3]. Rock microstructure intrinsically determines the rock storage capacity and flow capability [4][5][6]. Deformation of the rock microstructure influences the petrophysical properties of rock [7,8]. Accordingly, the relationship between permeability reduction and microstructure deformation should be investigated to understand the stress sensitivity mechanism [9,10]. Considering the nanoscale pore and fracture in tight sandstone and shale formations, a high-resolution microscale approach is required to further investigation of this.
Attempting microscale investigations, microscopic imaging technology combined with digital image correlation (DIC) techniques have been employed on various rock types. The microscale deformation and damaging experiments were usually conducted for investigating the fracture initiation and propagation mechanism. At the millimeter scale, the sandstone indentation experiments were carried out until failure, and the fracture length evolution was characterized by analyzing DIC results [11]. The experiment results showed that the DIC was efficient in tracing the deformation behavior and identifying the length of cracks formed and the onset of the tensile fracture. By implementing disk-shape experiments and the DIC method, He and Hayatdavoudi [12] identified that the direction of fracture initiation in sandstone grains and pore-filling materials is along the diametric loading plane. However, the above experiments have not reached the pore and fracture scales, which is expected to be critical for studying the stress sensitivity effect on the microlevel. For investigating the mechanical properties of mudstone and shale, Wang et al. [13][14][15] presented an innovative microscale experimental approach, which consists of in situ mechanical loading test in the environmental scanning electron microscopy (ESEM) and full field measurement by DIC techniques. The high resolution of the experiments was guaranteed by the use of ESEM. Our research work applied Wang's microscale experimental approach for studying the 2D microstructural deformation of tight sandstone but used SEM instead of ESEM. Most of the microscale deformation experimental studies served for studying the fracture initiation and propagation mechanisms or the mechanical properties of some certain rock types. The microscale deformation results in our work were innovatively combined with the digital rock technique for studying the microscopic stress sensitivity effects.
The development of digital rock physics (DRP) technology also provides opportunities for the microscopic stress sensitivity studies. On one hand, DRP technology has become one of the most powerful tools to calculate permeability by microflow simulation, and feasible calculation methods focusing on several main reservoir rock types have been established and recognized by industry [16][17][18]. For medium-high permeable sandstone and tight sandstone rocks, Yang et al. [19] extracted the pore network from digital rocks under different states of stress and calculated the corresponding permeability values by microflow simulation. The results showed that tight sandstones have more severe stress sensitivity than medium-high permeable sandstone. On the other hand, three-dimensional in situ micro-CT experiments have been used for the microscopic stress sensitivity studies. For example, such experiments for coal cores were carried out in Zhang's research work, and several 3D digital cores at different stress states were obtained. The cleat change law and the relation between cleat deformation and the coal permeability decrease were analyzed [20]. However, 3D in situ scanning and loading experiments are very expensive and time consuming. In addition, the integration of 3D scanning equipment and loading equipment can provide only a rela-tively low resolution of digital cores, down to several micrometers. For unconventional reservoirs such as shale and tight sandstone, pore diameters are at the nanometer scale [19][20][21][22][23]. To overcome the resolution limits and make the study time and cost effective, we proposed a morphological image processing method to simulate the deformation of digital cores with increasing stresses. A series of simulated cores at different stress states were used for pore network model extraction, and the corresponding permeability was calculated.
In this paper, we presented a novel workflow for studying the microscopic stress sensitivity of tight sandstone. To implement the workflow, the high-resolution 2D in situ SEM loading scanning experiments were firstly carried out to quantify the deformation at the pore and fracture scales. Then, the digital core deformation simulation was conducted and followed by microflow simulation to calculate the permeability changes under different stress states. The feasibility of the combination of in situ SEM studies and digital core deformation simulation was validated with the porositypermeability integrated measurements under triaxial stresses.

Materials and Methods
2.1. 2D In Situ Loading Scanning Experiment. To study how the micro/nanoscale pores and throats deform with loading, high-resolution 2D in situ loading scanning experiments were conducted, and high-resolution SEM images at different stress states were obtained. DIC was applied for image analysis; thus, the strains of the pore and throat areas were calculated.
The 2D in situ loading scanning experiment combines scanning and loading equipment. In this study, a Zeiss Sigma 500 field emission scanning electron microscope (FE-SEM) and Gatan Microtest 5000 loading equipment were used ( Figure 1). The accelerating voltage of this FE-SEM is 0.02 kV to 30 kV, and the resolution can be 0.8 nm at 15 kV and 1.6 nm at 1 kV; the magnification range is 12 to 1,000,000. The accelerating voltage of 8 kV was chosen for imaging acquisition in this study. The Zeiss Sigma 500 is equipped with an Inlens probe (secondary electron detector in the lens tube), an SE2 probe (secondary electron detector in the sample chamber), and a backscatter probe

Geofluids
(backscattered electron detector). The yield of secondary electrons depends mainly on the surface morphology of the sample. The yield should be high in regard to the convex part, and there should be no secondary electron generation in the pore and throat areas. The gray value of the image is the reflection of the yield of secondary electrons: if the yield is high, then the image gray value is high, causing the image to be bright; in contrast, the pores and throat area lead to a low gray value of the SEM image, so the image is black. To better observe the pore and throat structure of the tight sandstone, the image acquisition in this paper was completed with the SE2 probe. The Gatan Microtest 5000 in situ loading equipment can achieve uniaxial tensile and compression with a maximum force of 2000 N in the tensile mode and a maximum force of 5000 N in the compression mode. The loading velocity can be adjusted within 0.1~1.25 mm/min, and the data acquisition time interval is 100 ms~5 s. The force sensor accuracy is 1%, the displacement sensor accuracy is 10 μm, and the maximum sample size is 9 mm × 5 mm × 5 mm.
The tight sandstone sample used for the experiment is 9 mm in length, 5 mm in width, and 5 mm in thickness. The sample was collected from downholes at a burial depth of 3500 m in the Chunxiao Oil Field in China. The sample was cut into a suitable size by wire-electrode cutting to ensure the flatness and smoothness of the two lateral sides in contact with the loading equipment. To observe the pores and throats better and quantitatively characterize nano-micropores and throat structures, the flatness of the observation side must be guaranteed [24,25]. Consequently, the observation side was polished by the mechanical polishing method followed by the iron beam milling polishing method. After polishing, a carbon coating was applied to enhance the electrical conductivity of the sample surface, which is necessary for highresolution SEM image acquisition. After preparation, the rock sample was placed into the compression clamp of the in situ loading equipment, and then, the loading equipment was installed into the sample chamber of the FE-SEM, which was kept at high vacuum during the experiment. The loading equipment and the scanning equipment were controlled by two different computers.
The sample was loaded with uniaxial compression loading along the sample length direction. The loading rate was 0.1 mm/min, which is the minimum loading rate of the facility. In this in situ loading scanning experiment, images were acquired at 0 MPa, 6.5 MPa, and 10.5 MPa. To avoid rock debris generated during rock breakage, this experiment was not loaded to higher stress states. When acquiring images every time, loading was suspended, and the loading equipment was in measurement mode in which the load was constant. Since the image acquisition lasted for several minutes, maintaining a higher loading would have increased the image noise.
The SEM images acquired by the 2D in situ loading scanning experiment can be analyzed by the DIC method, and then, the specific displacement and strain information of the microstructure can be obtained. The DIC method is an optical method that calculates the displacement and strain of the points all over the image by comparing the initial and deformed images [15,26]. The working theory of DIC is briefly introduced here. First, a subset is chosen from the initial image, and then, the location of the subset in the deformed image is searched by correlation coefficient optimization. The correlation coefficient is defined by the grayscale gradient of the subset, and the optimal value of the correlation coefficient can indicate the location of the subset in the deformed image. Then, the calculation of the displacement and strain fields can be accomplished by comparing the subsets in the initial and deformed images [27]. In this study, open source DIC software was used for SEM image analysis   3 Geofluids and calculation of strain maps [28]. Compressive strain was considered negative, and extensional strain was considered positive. The Green-Lagrangian strain was used for calculation, which is defined by the following equation: where E xx means the horizontal strain, u means the horizontal displacement, and v means the vertical displacement.
2.2. 3D Digital Core Imaging Experiment. Micro-CT imaging experiments are needed for high-resolution 3D digital core acquisition, which is the foundation of 3D digital core deformation simulations. The same sample as in the 2D in situ loading scanning experiment (with a diameter of 2 mm and length of 3 mm and sampled by a special drill bit) was used for the 3D digital core imaging experiment. With the Xradia CT-520 Versa from Zeiss (Figure 2), 1003 images with 992 × 1024 pixels were obtained. The resolution of the images is (0.7 μm) 3 for each voxel. Considering the calculation costs, 400 images were selected and then cropped to 400 × 400 × 400 pixels (Figure 3) for the digital core deformation      5 Geofluids simulation. Images were filtered by the median filter method to reduce image noise [29]. The CT images are gray images, which need to be segmented. Segmentation is a vital step for the following analysis; here, the Otsu method was used for segmentation [30].

Porosity-Permeability Integrated Measurements.
Porosity-permeability integrated measurements were performed to obtain the rock porosity and permeability at different stress states under real compression conditions. The porosity values were used as the constraints of the 3D digital core deformation simulation, and the permeability values were regarded as the reference data for evaluating the simulation reliability. The porosity-permeability integrated measurements were carried out with samples 25 mm in diameter and 20 mm in length, obtained from the same source as in the former experiments. A Poro-Perm PDP-200 Pulse Per-meameter was used for measurements. Porosity was measured with helium, and permeability was measured with nitrogen. The permeability was tested with the pulse decay method.
A series of values of porosity and permeability were obtained with increasing confining pressures from 5.5 MPa, and the pore pressure was constant at 5.5 MPa. Effective stresses were calculated with the following equation according to Terzaghi effective stress theory [31]: where P eff is the effective pressure, P c is the confining pressure, and P p is the pore pressure. In this paper, five stress states of 0 MP, 6.   Geofluids and permeabilities are shown in Figure 4. The permeability decreased sharply to half of the original value, especially at the former two stages, while the porosity slightly decreased with increasing effective stress.
2.4. Simulation of 3D Digital Core Deformation. 3D rock in situ loading scanning experiments are expensive and usually have low resolution. To efficiently obtain high-resolution digital cores at different stress states, the 3D digital core obtained from the micro-CT imaging experiment and porosity values from porosity-permeability integrated measurements were combined to simulate the microstructure deformation of the 3D digital rock. Digital cores at different stress states could be acquired by simulation, and then, pore network models were extracted, with which permeabilities and pore throat parameters could be calculated.
We proposed an image processing method to simulate the microstructure deformation of digital rock. When rock is under confining pressure, the rock particle edge extends forward, and the pore edge shrinks toward the inside. As Figure 5 shows, the radius of spherical particles increased, and the edge extended forward with increasing loading stress [32]. For the rock under uniaxial loading, the pore edge shrank along the loading direction. Based on the fundamental theories of morphology, the algorithm of image erosion can follow the loading. Our method used related MATLAB functions to process the images of digital rocks and to simulate the pore structure changes ( Figure 6) under different states of stress. Porosities under different stresses tested by porosity-permeability integrated measurements can be used as the constraint parameters for such simulations.
Images obtained from the 2D in situ loading scanning experiment were used for illustration. The image at 0 MPa   9 Geofluids was segmented into a binary image as the initial image for deformation simulation. The porosity at 6.5 MPa was the simulation constraint. The erosion simulation was accomplished by MATLAB, and then, the deformation image was obtained (Figure 7), where black corresponds to pores and throats and white corresponds to the rock skeleton. When comparing the resulting three images, the real deformation image shows that stress closes some pores (Figure 7), and the simulated image shows an identical scenario to the real case.
The pore diameter and pore shape factor are calculated from Figures 7(b) and 7(c) and are shown in Figures 8 and  9. Here, the 2D pore shape factor is calculated by the following equation.
where L P is the pore perimeter and A p is the pore area. A comparison of the latter figures indicates that the statistical parameters of the real deformation image and simulation image exhibit high similarity. The range of pore diameters is the same, and the frequency distributions both indicate that small pores are the dominant component of pores and that their frequency decreases with increasing effective stress. The range of the pore shape factor is similar, and its frequency distribution demonstrates that most pores are elliptical, while some exceptionally large shape factors indicate the presence of long and narrow throats. According to all these features, it can be concluded that the simulation method reproduced the real deformation image very well.
In this study, we combined the 3D digital core obtained from the CT experiment and the porosities measured from the porosity-permeability integrated measurements for the simulation. The porosities were considered as the simulation constraints. The parameters of the erosion algorithm were adjusted several times until the porosity of the eroded digital core was equal to the measured porosity. By repeating the above steps, a series of simulated deformed digital cores at each stress state were obtained.
Pore network modeling, the typical method for pore and throat characteristic parameter analysis, is also used for permeability calculations [33,34]. In the present work, pore network models are extracted from the simulated deformation digital rocks with a modified maximal ball algorithm [35]. Pore and throat parameters of simulated deformation digital rocks were calculated from pore network models. Singlephase flow simulations using pore network models were carried out to obtain the absolute permeabilities, with the hypothesis that the flow obeys Darcy's law, and the simulation was implemented by our in-house program.

Results of the 2D In Situ Loading Scanning Experiment.
The SEM images at different stress states from the 2D in situ loading scanning experiment were investigated with the DIC method to calculate the strain maps. In this study, the loading direction was taken as the x-axis, as the coordinates in Figure 10(a) show. Considering the uniaxial loading scenario, only E xx strain was demonstrated. The original SEM images and contour strain maps were overlapped to clearly view the strain at different positions. One of the typical observation zones was selected for strain study during the loading process ( Figure 10). Figure 10(a) shows that there are many pores and throats in the observation zones. The strain maps indicate a highly heterogeneous strain distribution over the image area. As shown in Figure 10(b), the throats perpendicular to the loading direction were compressed and caused larger strains than the pores. The strain contour maps also demonstrated the strain concentration around the throat areas, where the strain contour lines were denser. The strain concentration became more severe as the effective stress increased (Figure 10(c)).

Geofluids
The study of strain maps during loading indicated that the deformation of throats was more intense than that of pores.

3.2.
Results of 3D Digital Core Deformation Simulation. The digital core obtained from the micro-CT imaging experiment with 400 × 400 × 400 voxels was selected for deformation simulation. Note that the calculated porosity of the original digital core by Avizo is 17.2%, which is different from the tested porosity of 7.12% at 0 MPa by porosity-permeability integrated measurement. The digital core is a 0:28 × 0:28 × 0:28 mm 3 cube, while the rock sample used in the integrated measurement is a cylinder of 25 mm diameter and 20 mm length. Considering the large size differences of the digital core and the laboratory sample, the different porosity values should be attributed to the scale effect [36,37]. Therefore, the porosity change rate rather than porosity in porositypermeability integrated measurements was regarded as the simulation constraint. With the deformation simulation method introduced before, a series of simulated deformed digital cores were obtained at different stress states ( Figure 11). The white part represents the pores and throats, and the red part is the rock skeleton. Notably, the white part slightly changes with increasing stress. By extracting pore networks from the simulated digital cores, the absolute permeability of each digital core was calculated, as shown in Figure 12. The laboratory-tested permeabilities are different from the simulated values, but the changing trends are very similar. Accordingly, the reliability of the deformation simulation method was demonstrated again.
The radii of the pores and throats were calculated from the extracted pore network models, and the radius frequency distributions of the pores and throats are presented in Figures 13 and 14, respectively. From Figure 13, we found that the total pore number decreased from 8752 at 0 MPa to 5894 at 21 MPa. The pore radius frequency distributions also varied with increasing effective stress. From 0 MPa to 6.5 MPa, the frequency count of pores smaller than 1 μm decreased, while the frequency of larger pores was relatively higher, which indicated that the small pores diminished first at the beginning of loading. From 6.5 MPa to 21 MPa, the pores smaller than 1 μm showed increasing frequency, while the frequency of the larger pores decreased. The analysis revealed that the larger pores diminished into small pores more intensely during the latter period of loading. In addition, the range of pore radius also varied with increasing stress. As the effective stress increased from 0 MPa to 21 MPa, the maximum pore radius shrank from 17 μm to 15.1 μm, which demonstrated the shrinkage of the largest pores with increasing effective stress. Figure 14 shows that the total number of throats decreased even more severely with increasing stress, from 16634 at 0 MPa to 8573 at 21 MPa. The throat radius frequency distributions varied with increasing stress. From 0 MPa to 6.5 MPa, the frequency of throats whose radii were smaller than 0.8 μm decreased, the frequency of throats whose radii were larger than 11 μm dropped to zero, and the frequency of throats whose radii were between 0.8 μm and 11 μm increased relatively. Such variations demonstrated that the largest throats diminished into smaller throats during the early loading period. From 6.5 MPa to 21 MPa, the frequency of throats whose radii were smaller than 0.8 μm continued to increase, while the frequency of larger throats decreased, especially the largest throats. This is because increasingly larger throats were diminishing into smaller      throats. In addition, the range of the throat radius changed more severely than the pore radius. The maximum throat radius shrank from 12 μm at 0 MPa to 8 μm at 21 MPa. Based on the above results and analysis, it can be concluded that throats shrink more intensely than pores as the effective stress increases.

Conclusion
In this paper, a 2D in situ loading scanning experiment was carried out to investigate the microstructure deformation of tight sandstone rock samples. The DIC method was applied to SEM images to calculate the strain maps during uniaxial loading. The results showed that throats perpendicular to the loading direction exhibited strain concentration, and the strain in those areas was obviously higher than that in the pore area.
Additionally, we proposed a deformation simulation method based on morphological image processing. The statistical parameters between real deformation images and simulation deformation images were analyzed to validate the simulation deformation method.
With the simulation deformation method, a digital core obtained from the micro-CT imaging experiment was used for the 3D digital core deformation simulation. The porosities tested from porosity-permeability integrated measurements were used for simulation constraints. Pore networks were extracted from the simulated deformation digital cores, and then, the permeabilities could be calculated by pore network modeling. The simulated permeability showed a similar changing trend as the laboratory-tested permeability. During the loading process, the frequency counts and range of pore and throat radii were presented and analyzed. The results provide strong evidence that throat deformation was more severe than pore deformation under stress, which is considered the main reason for the severe stress sensitivity effect in tight sandstone reservoirs.

Data Availability
The experiment data used to support the findings of this study are currently under embargo while the research findings are commercialized. Requests for data, 12 months after publication of this article, will be considered by the corresponding author.

Conflicts of Interest
The authors declare that they have no conflicts of interest.

Acknowledgments
This study was supported by the National Natural Science Foundation of China (51474224).    16 Geofluids