Porosity and Permeability Evolution with Deviatoric Stress of Reservoir Sandstone: Insights from Triaxial Compression Tests and In Situ Compression CT

Shandong Provincial Key Laboratory of Marine Environment and Geological Engineering, Ocean University of China, Qingdao, China Laboratory for Marine Geology, Qingdao National Laboratory for Marine Science and Technology, Qingdao, China Key Laboratory of Marine Environment and Ecology, Ministry of Education, Ocean University of China, Qingdao, China LML CNRS FRE 3723 and Centrale Lille, Villeneuve D’Ascq Cedex, France


Introduction
Underground natural gas storage projects are developed into three different geological structures: depleted gas or oil reservoirs, aquifers, and salt caverns [1][2][3]. For storage projects in depleted sandstone reservoirs, porosity and permeability are the two most important characteristics [4]. Porosity will influence the capacity to store natural gas, and permeability will influence its deliverability rate.
Under in situ conditions, stress variation during gas injection into reservoir rocks will cause rock deformation and change the porosity and permeability. Better knowledge of porosity and permeability properties of geomaterials requires understanding of the coupling between deformation and pore fluids [5][6][7][8][9][10][11]. For high-porosity (>15%) reservoir sandstones, the failure mode will undergo a transition from brittle faulting (at low confining stresses) to ductile flow (at high confining stresses) under triaxial compression conditions [12]. Permeability decreases continuously at both of the two different failure modes due to compression of microcracks and increase in tortuosity [13][14][15][16][17]. For brittle faulting, porosity firstly decreases in the volumetric compression stage and then increases during volumetric dilation, while the porosity will continuously decrease for ductile flow [18][19][20][21].
The variation of porosity and permeability under hydrostatic or deviatoric loading conditions is thought to be related to the variation of microstructures. Digital rock physics has been used in the past decades to understand pore-scale processes and estimate macroscopic rock properties [22]. Computed tomography (CT) is a powerful nondestructive scanning technique. It can provide important insights to characterize the internal nano/micro/mesoscale structures of many geomaterials [23][24][25][26][27]. In situ compression CT, which scans internal structures of materials under progressive loadings, is now feasible. The structural damage and fracture evolution can be examined and related to the loading process [25,[28][29][30][31][32][33][34]. Fonseca et al. investigated the micromechanisms of inelastic deformation in Fontainebleau sandstone by a triaxial in situ compression CT [28]. Renard et al. quantitatively revealed the high-resolution 3D information about damage evolution of marbles undergoing brittle failure with an in situ X-ray transparent triaxial deformation apparatus [35]. McBeck et al. captured the path to macroscopic shear failure within shale using time-resolved in situ X-ray CT [36]. Understanding material deformation and damage on the sub-micrometer scale is useful to explain the evolution of porosity and permeability [37]. Based on CT scanning, some pore-scale modeling methods can be used to simulate and explain transport behavior of geomaterials [38][39][40]. The image-based pore-scale modeling methods are mainly categorized into two classes: topology-central methods and morphology-central methods [41]. Based on the extracted pore network, a series of transport simulation methods are developed such as quasistatic fluid modeling, dynamic flow modeling, solute/colloid transport, and reactive flow modeling [42].
In this study, the porosity and permeability evolution with deviatoric stress of reservoir sandstone was investigated by triaxial compression tests and a uniaxial in situ compression CT test. Common methods to investigate permeability evolution with deviatoric stress need to apply a confining stress, such as 5 MPa, which will lead to the compression of microstructures. The applied pressure is to ensure that pore fluid passes through the sample. With our experimental technique, we can measure gas permeability at a very small confining stress (0.2 MPa), which is very close to uniaxial compression conditions. The interest in permeability under low confinement is because the CT test was performed under uniaxial conditions. The uniaxial in situ compression CT test was made in order to better explain the evolution of porosity and permeability with loading from the point of view of the microscale, while the trade-off between sample size and imaging resolution still limits the detection and visualization of induced microstructural changes [43]. Therefore, we chose a small sample size for the uniaxial in situ compression CT test. The pore properties were analyzed with Avizo 9.0. A Maximal Ball Method (MB), which is a morphology-central method, was then used to extract the pore network based on CT images, and the absolute permeability evolution with deviatoric stress was simulated. The failure mode of the sandstone in this experimental study is brittle faulting. Porosity measurement needs to inject gas in the sample. It is impossible to measure porosity at very low confining pressures. Therefore, porosity evolution with deviatoric stress was measured at different confining pressures from 5.5 to 30.5 MPa.

Test Material and Apparatus
2.1. Test Material. The reservoir sandstone used in this research is obtained from the Vosges Mountains in the east of France. The sample size for the triaxial compression test is 70 mm in height and 37 mm in diameter. The size of the sample for the in situ compression CT test is 12.14 mm in height and 5.31 mm in diameter. All samples were in dry conditions, and the tests were conducted at a room temperature. The grain size is mainly in the range between 250 and 300 μm. Figure 1(a) is the CT image of the sample for the triaxial compression test (70 mm in height and 37 mm in diameter). We can observe several small sedimentary bands. The approximate position to drill the small sample for the uniaxial in situ compression CT test is also shown in Figure 1(a). Figure 1(b) is the CT image of the in situ compression CT. Figure 1(c) is the pore structure of the same vertical cross section extracted with Avizo 9.0. It is observed that the pore structure is uniform and most part of the small sample does not pass the sedimentary bands.

2.2.
In Situ Compression CT. The uniaxial in situ CT data was cited from the work of Hu et al. (2018). The axial stress was applied by a press machine. Four different scans were made at different deviatoric stresses (0. 45, 22.50, 31.50, and 36.00 MPa). During scanning, the deviatoric stress was kept constant. The loading machine rotated around the rotating stage. After one scan was finished, deviatoric stress was increased to another level. The deviatoric stress loading speed applied by the loading machine was 0.12 mm/min. The initial small stress (0.45 MPa) was to prevent the sample from moving during machine rotation. At each loading step, the total scan time was about 80 min. The acceleration voltage was 110 kV. The current was 38 μA. The resolution was 5 μm. Each scan generated 2742 projections.

Triaxial Compression
Apparatus for the Porosity and Permeability Test. Figure 2(a) is a servo control triaxial compression apparatus. Two strain gauges were glued to the sample at midheight with an angle of 180°to each other. The gauge locations were first applied a thin layer of fast-cured epoxy resin to seal the surface pores [44]. Each gauge can measure one axial (ε a ) and one lateral (ε l ) strain at the same point. The final axial and lateral strains are the average of the two gauges. Volumetric strain (ε v ) is ε a + 2ε l . The surface of the sample was sealed with a membrane to avoid the leak of confining fluid into the sample. The strain values were recorded using a program code in LabVIEW. A loading machine (Zwick Roell) was used to apply the deviatoric stress (σ 1 -σ 3 ). The loading speed was 0.12 mm/min. Pure argon was used as the pore fluid in our test. Porosity variation with deviatoric stress was measured at four different confining stresses (5.5, 10.5, 20.5, and 30.5 MPa). The procedure measuring porosity variation with deviatoric stress is as follows: (1) All valves were firstly closed. Then, about 0.8 MPa (P 1 ) gas pressure was injected into the accumulator (the blue part in Figure 3(a)). The volume of the accumulator (V 1 ) is 46.317 ml.
(2) The designed confining stress was applied. When P 1 became stable, valves 1 and 3 were opened. The current gas pressure became P 2 . The volume of the green 2 Geofluids part (V 2 ) is 6.864 ml. When P 2 became stable, deviatoric stress was applied. The manometer is connected to a computer, which can record the variation of gas pressure every second.
(3) The porosity at different times can be calculated with the ideal gas law: P 1 V 1 = P 2 ðV 1 + V 2 + V pore Þ. Therefore, the volume of pores (V pore ) at different deviatoric stresses can be calculated. Porosity is the ratio between the pore volume and the sample volume.
A permeability test was made at different deviatoric stresses at a confining stress of 0.2 MPa. The stress state is very close to a uniaxial compression condition. When measuring permeability, valves 1, 2, and 4 were opened and valves 3 and 5 were closed. The permeability of this sandstone is very high, and the pressure loss at the downstream cannot be ignored. The best technique for measuring permeability is to maintain a constant flow rate while measuring the gas pressure [44]. Therefore, the cell must be accurately calibrated prior to fluid injection. This calibration was made in the absence of any sample inside the cell, as shown in Figure 2(b). The flow rate was set as 0.15 Nl/min. The loss of pressure had an average value (measuring 11 times) of 0.0746 MPa. After deviatoric stress was loaded to a constant value, gas was injected into the sample at a constant flow rate of 0.15 Nl/min. About 5 minutes later, the injection pressure was stable and the measured pressure was recorded. The following equation based on Darcy's law is used to calculate the permeability: where K is the permeability (m 2 ). L is the height (length) of the sample (m), q is the flow rate (Nl/min is the case for our flowmeter), μ is the viscosity of argon (2:2 × 10 −5 Pa · s). S is the cross-sectional area of the sample (m 2 ). P a is the atmosphere pressure (P a ), P i is the measured injection pressure, and P 0 is the calibrated gas pressure (P a ).

Imaging Process
The scanned images were processed and analyzed by Avizo 9.0. This software is a 3D visualization, analysis, and modeling system [45]. The processing procedure is as follows: (1) Because the sample was compressed and rotated during scanning, the images generated from CT scanning at different loading stages will change and it was impossible to find and compare the same crosssectional area at different loading stages. The process of "register" was first applied to process the images. The images of the second, third, and fourth scan were registered with the first scan. This process will choose 3 Geofluids the middle horizontal and vertical cross section of the first scan as the reference. After "register," we can get one same horizontal cross section and two same vertical cross sections in the middle position of the sample at four different loading stages.
(2) A median filter was applied to reduce the image noise.
(3) A grayscale threshold was used to segment the solids (77-255) and pores (0-77). The same threshold was chosen for all the different scans.
(4) There will be plenty of noisy points at the two ends of the sample touching with the loading plates. Therefore, a subvolume (cuboid) was extracted to delete the two parts near the loading plates.
(5) An envelope and erosion technique was then applied to obtain a cylindrical structure of mask (total volume of the digital sample). Then, the final cylindrical pore structure can be obtained by the intersection between the cylindrical mask and the cuboid pore structure obtained in step (3).

Pore Network Extraction and Absolute Permeability Simulation
The permeability simulator needs a cubic model. So, we extracted a subvolume from the cylindrical sample. Figure 4 shows the three middle reference cross sections and the extracted subvolume for pore network modeling.

Geofluids
The extracted subvolume can be transformed into a 3D raw format in Avizo 9.0. The Maximal Ball Method (MB) was used to extract the pore network [46]. In this method, a MB touches grain surfaces and is not a subset of any another MB. The maximum balls are defined as pores. The chains of balls between pores are defined as throats. A two-step searching algorithm was used to search the nearest solid to define a void ball. A clustering process to define pores and throats was invented by affiliating the maximal balls into family trees according to their size and rank [46]. The porosity of the pore networks extracted by the MB method is 16

Geofluids
shows the pore network extracted by the MB method. This pore network includes pores and throats.
A pore-scale simulator developed by Valvatne was used to predict the absolute permeability of the pore networks [47]. Fluid is assumed to be incompressible. The absolute permeability of the pore network is derived from Darcy's law: where μ is the viscosity of a single-phase fluid (cp). q tsp is the total single-phase flow rate (m/s). L is the length of the pore network (m). S is the cross-sectional area (m 2 ), and P inlet − P outlet is the pressure drop (P a ). The total single-phase flow rate is found by solving for the pressure everywhere and imposing mass conservation at every pore i.
where j runs over the whole throats connected to pore i. Fluid flow is assumed to be incompressible, and pressure drops to flow are insignificant compared to the capillary pressure.
The flow rate q p between two pores i and j is as follows: where g p is the fluid conductance. l is the length between the pore centers. Ф p is the phase potential.

Experimental Porosity Evolution with Deviatoric Stress.
When valves 1 and 3 were opened, the pressure (P 2 ) was around 0.55 MPa. Based on Terzaghi's effective stress principle, the effective confining stress is around 5, 10, 20, and 30 MPa. The stress-strain relationship at different confining stresses is shown in Figure 5. Therefore, the occurrence of volumetric dilation will be delayed with the increase in confining stress.
The unloading modulus is shown in Table 1. The modulus increases gradually with the increase in deviatoric stress. When approaching the peak stress, the increase rate becomes very small or the modulus even decreases in some cases. Therefore, we can emphasize that microcracks will be compressed in the axial direction. At low confining stress, the increase in unloading modulus is more obvious than that at high confining stress. The reason is that initial pores and microcracks are closed due to high confining stress. Figure 6 shows the relationship between porosity and volumetric strain. We can find that the relationship is nearly linear. At confining stresses of 5.5 and 10.5 MPa, the initial porosity (φ 0 ) is smaller than the last measured porosity (φ l ) before failure. At a confining stress of 20.5 MPa, φ 0 is nearly equal to φ l . And at 30.5 MPa, φ 0 is much bigger than φ l . Therefore, the level of volumetric dilation decreases with the increase in confining stress.
The porosity we used is the Lagrangian porosity which refers the current porous volume to the initial sample volume [48]. If the solid matrix is assumed to be nondeformable, the deformation mainly comes from the change in the pore space. Therefore, the porosity can be expressed as [18] where φ is the current porosity, φ 0 is the initial porosity, ε v is the current volumetric strain, and e 0 is the initial void ratio. The initial porosities we measured were 18.23%, 18.81%, 18.43%, and 18.28%, respectively, at confining pressures of 5.5, 10.5, 20.5, and 30.5 MPa. Figure 7 shows the measured and calculated porosity evolution with deviatoric stress. It is obvious that the porosity calculated with equation (4) is in good agreement with the measured results. At the compression stage, porosity decreases with the increase in deviatoric stress and increases at the dilation stage. Although the measured and calculated results are in good agreement, it should be noted that the porosity variation is very small. The last porosity (φ l ) values we measured before failure are 18.39%, 18.90%, 18.43%, and 18.24%, respectively, at confining pressures of 5.5, 10.5, 20.5, and 30.5 MPa. After failure, the residual deviatoric  .46%, respectively. The porosity will increase (φ r − φ l = 1:4%, 0.59%, 0.38%, and 0.22%) after failure due to the development of macrocracks. The porosity increase (φ r − φ l ) after failure decreases gradually with the increase in confining stress. We can conclude that the volume of induced cracks by deviatoric stress at low confining stress is bigger than that at high confining stress. If the compressibility of the solid matrix is taken into account, the variation of the Lagrangian porosity can be deduced as follows [48]: where φ is the current porosity, φ 0 is the initial porosity, and ε v is the volumetric strain of the skeleton (the variation of volumetric strain is positive when the sample is compressed). ε s is the volumetric strain of the solid matrix: where σ s is the mean stress applied on the solid matrix (MPa) and K s is the solid matrix bulk modulus (GPa). The values of K s we measured at confining stresses of 5, 10, 20, and 30 MPa are 35, 39, 41, and 44 GPa, respectively [49]. The value of σ s where σ m is the mean stress applied on the skeleton (MPa). If K s is assumed to be constant during applying deviatoric stress in the compression stage, we can estimate the porosity evolution with deviatoric stress. Figure 8 shows the calculated porosity with equation (6) and the measured porosity in the compression stage. Trapped porosity can be estimated by the difference between the calculated and measured results. For this sandstone, the trapped porosity is very small, which can be neglected during triaxial compression. Therefore, the pore structure of this sandstone is nearly all connected. Increase in deviatoric stress or confining stress will only compress or close the pore structure but not induce occluded pores. The response of pore structure to ambient stress of this high-porosity sandstone is quite different from that of tight sandstone with a porosity of about 5% [50]. For this tight sandstone, some big pores will be trapped inside the solid matrix with the increase in confining stress from 3 to 20 MPa.

Porosity Evolution with Deviatoric Stress Based on Image
Analysis. After the images are processed by the method mentioned in Section 3, some useful pore properties can be obtained. "Label analysis" in Avizo 9.0 is used to analyze the pore properties, which can recognize different pores and calculate the volume of these pores.
If we obtain the volume of pores and mask, the Eulerian porosity, which refers the current pore volume to the current sample volume [48], can be calculated. The pore properties are listed in Table 2. When deviatoric stress increases to 22.50 MPa, porosity decreases from 16.19% to 16.08%. With continuous increment in deviatoric stress, porosity increases to 16.12% and 16.18%, respectively. The volumetric compression and dilation are successfully captured by the in situ compression CT test. We can predict that 22.50 MPa is around the transition from volumetric compression to dilation. For this sandstone, under uniaxial compression condition, it is normal that the variation of volumetric strain is not as significant as that at high confining stresses [51]. The maximum volumetric strain of this sandstone is only 0.000603 under uniaxial compression. And the small variation of porosity with deviatoric stress is consistent with volumetric strain variation under the uniaxial compression test.   Figure 6: Relationship between measured porosity and volumetric strain.

Geofluids
The variation of trapped porosity with deviatoric stress obtained by using a digital rock sample is also shown in Table 2 and Figure 9. The trapped porosity is very small. When deviatoric stress increases from 0.45 to 22.50 MPa, the decrease is very small from 0.122% to 0.116%. The reason may be the closure of some initial pores. At the volumetric dilation stage, it is nearly constant at 0.115%. These results are also consistent with that estimated by triaxial compression tests. The experimental technique used in this study makes it possible to determine permeability at very low confining stresses. The Reynold number we calculated based on equation (9) is around 4 under our test conditions, which indicates that the flow condition is laminar flow. Table 3 shows the measured relative injection pressure at different deviatoric stresses. The maximum injection gas pressure is smaller than half of the confining pressure (0.2 MPa). The difference between confining pressure and injection gas pressure will keep the measured permeability results reliable.

Permeability Evolution with Deviatoric Stress
where Re is the Reynold number. d is the hydraulic diameter (m). q is the flow rate (m 3 /s). κ is the kinematic viscosity (m 2 /s). S is the cross-sectional area of the sample (m 2 ). Permeability evolution with deviatoric stress is represented in Figure 10. When deviatoric stress increases from 0 to 10.14 MPa, the decrease in permeability is very significant: from 6:18 × 10 −13 to 1:22 × 10 −13 m 2 . The significant decrease is due to the closure of initial pores and microcracks. With the increase in deviatoric stress, permeability decreases continuously because of the compression of pores and microcracks. At deviatoric stresses of 30.36 and 40.01 MPa, the sample is already in the volumetric dilation (observation from the relationship between volumetric strain and deviatoric stress). However, permeability still decreases. For highporosity sandstone (>15%), the decrease in permeability at   9 Geofluids the volumetric dilation stage is due to the increase in tortuosity [14][15][16]52]. Sand particles will be produced at this stage. Figure 11 shows the photograph of the failed sample. Sand production phenomenon and microcrack evolution are also directly observed by our in situ compression CT test as shown in Figure 12. The size of each figure is 2:5 × 2:5 μm 2 . a1-h to d1-h and a2-h to d2-h are the same sand particles in the middle horizontal cross section. a3-v to d3-v and a4v to d4-v are the same sand particles in the middle vertical cross section at different deviatoric stresses. Microcracks are indicated by dotted lines. The microcracks between different sand particles in the horizontal direction will be compressed.
At a deviatoric stress of 36 MPa, some microcracks in the vertical direction inside some sand particles have developed. These microcracks (or smaller sand particles) will increase the tortuosity and lead to the decrease in permeability at the volumetric dilation stage. Equation (10) is the Katz-Thompson equation. This equation considers the effect of tortuosity on permeability [53]. During volumetric dilation, both porosity and tortuosity increase. From our experimental results, we can conclude that the influence of tortuosity increase on permeability is more significant than that of porosity increase. The decrease rate of permeability (difference of permeability divided by difference of deviatoric 18 where τ is the average pore network tortuosity and d c is a throat diameter (m).
The permeability damage coefficient can also be used to define the sensitivity damage degree induced by deviatoric stress [54].
where D is the permeability damage coefficient. K D is the permeability at different deviatoric stresses. K 0 is the initial permeability.
The damage coefficient at deviatoric stresses of 10.14, 20.04, 30.36, and 40.01 MPa is 0.803, 0.840, 0.851, and 0.856, respectively. The coefficient increases with deviatoric stress. From 0 to 10.14 MPa, the permeability damage coefficient increases to 0.803. The influence of damage induced by   11 Geofluids deviatoric stress on permeability is very significant at this stage. After that, the increase in the permeability damage coefficient is comparatively moderate.

Simulated Absolute Permeability Based on Pore Network
Modeling. Figure 13 shows the simulated absolute permeability based on pore network modeling. Permeability decreases from 1:174 × 10 −12 to 1:086 × 10 −12 m 2 when deviatoric stress increases to 36.00 MPa. The permeability is one order of magnitude bigger than our macroscopic experiment results. The tendency of permeability variation with deviatoric stress is similar to that in the macroscopic test. The decrease rates of permeability from the initial state to 36 MPa are 0:032 × 10 −13 , 0:0089 × 10 −13 , and 0:022 × 10 −13 m 2 /MPa, respectively. The permeability decrease rate is similar to the macroscopic experimental results. The permeability damage coefficient is 0.060, 0.066, and 0.075, respectively, which is smaller than that in macroscopic tests.
The simulated absolute permeability is about one order of magnitude bigger than the measured results. The damage permeability coefficient based on pore network modeling is also smaller than the measured results. On the one hand, the CT resolution is around 5 μm and the microcracks smaller than 5 μm may not be distinguished. Some Existed sand particles after failure Figure 11: Photograph of the failed sample at a confining stress of 0.2 MPa. 12 Geofluids   Figure 13: Simulated permeability evolution with deviatoric stress.
microcracks induced by deviatoric stress may not be reflected in the extracted pore networks. The fluid migration path in the extracted pore network may not be as complex as that in the actual rocks. The induced microcracks will also increase the tortuosity and decrease the permeability. On the other hand, the sedimentary band mentioned in Section 2.1 will also have an influence on the permeability. The macroscopic sample is slightly heterogeneous, and there are several sedimentary bands, which will reduce the permeability. As shown in Figure 1(b), the small-cored sample for CT does not completely cross sedimentary bands.

Conclusions
The evolution of porosity and permeability with deviatoric stress of reservoir sandstone is investigated by traditional triaxial compression tests and an in situ compression CT test. Volumetric dilation will be delayed with the increase in confining stress. Porosity firstly decreases with the increase in deviatoric stress due to compression of microcracks and pores. Then, it will increase at the volumetric dilation stage. For this high-porosity sandstone, nearly all pores are connected and trapped porosity is very small. Image analysis based on in situ compression CT supplies another insight to explain the porosity evolution with deviatoric stress. Based on image analysis, porosity variation with deviatoric stress is consistent with triaxial compression tests, and the trapped porosity is very small. Permeability evolution has been measured at low confining stress (0.2 MPa), which is close to uniaxial compression conditions. At the initial loading stage, permeability decrease is very significant due to closure of initial microcracks. The closure of microcracks in the horizontal direction can be directly observed by the in situ compression CT test. At the volumetric dilation stage, the decrease in permeability is due to the sand production and the increase in tortuosity. Observation of in situ compression CT confirms that sand particles may be fractured and sand is produced, and microcracks may be induced inside sand particles in the vertical direction at high deviatoric stress. The decrease in absolute permeability with deviatoric stress is confirmed by estimations based on pore network modeling.
Abbreviations ε a : Axial strain ε l : Lateral strain ε v : Volumetric strain ε s : Volumetric strain of the solid matrix σ m : Mean stress σ s : Mean stress applied on the solid matrix σ 1 -σ 3 : Deviatoric stress τ: The average pore network tortuosity μ: Viscosity of argon φ 0 : Initial porosity φ: Current porosity κ: Kinematic viscosity φ r − φ l : Porosity increase Ф p : Phase potential d: Hydraulic diameter d c : Throat diameter D: Permeability damage coefficient e 0 : Initial void ratio g p : Fluid conductance K: Permeability K 0 : Initial permeability K D : Permeability at different deviatoric stresses K s : The solid matrix bulk modulus l: Length between the pore centers L: Height of the sample P a : Atmosphere pressure P i : Measured injection pressure P 0 : Calibrated gas pressure P inlet − P outlet : Pressure drop in permeability simulation q tsp : Total single-phase flow rate in permeability simulation q: Flow rate Re: Reynold number S: Cross-sectional area of the sample.

Data Availability
Data is available based on request.

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