Porosity and permeability are two of the most important rock properties parameters to characterize the ability to store and transport hydrocarbons in porous media^{[1,2,3,4,5,6,7]}. Porosity and permeability usually decrease with increasing stress in both consolidated and unconsolidated porous media^[8,9]. The trend of porosity and permeability changes during the production processes depends on chemical, physical, and mechanical factors^[10]. Therefore, since carbonate reservoirs are typically heterogeneous and anisotropic, the full understanding of the changes of these two parameters plays a significant role in determining strategies for hydrocarbon production^[11,12,13]. Karacan et al.^[3] investigated the porosity and structural changes in sandstone porous media under triaxial loading conditions to measure the relationship of porosity and permeability. They concluded that in the ductile failure condition, the porosity was reduced by applying pressure in all directions, and it increased after approaching to ductile-plastic deformation condition. There was also a gradual reduction in the permeability due to loading; with the onset of fracturing, an increase in permeability was observed. If the stress is applied repeatedly, which is also called ''cyclic loading'', the damage to the rock sample will increase gradually. The rock will deform correspondingly, and may experience fatigue failure and decrease in strength^[14,15]. Teklu et al.^[16] investigated the overburden permeability and hysteresis of cores from different carbonate reservoirs, and proposed cyclical changes of permeability during pressure loading and unloading cycles. They found that the permeability decreased with the increase of stress and the process had some hysteresis.

The development of imaging technology in recent decades, which is referred to as digital rock physics (DRP), improves the ability to evaluate the reservoir rock properties and examine the internal structure of the rock^[17,18,19]. Vinegar and Wellington^[20] used computed tomography (CT) to monitor the oil production process for the first time. They described the rock physics using CT scanning and X-ray methods and proposed a protocol for multiphase fluid flow. Afterward, many researchers used CT scanning to measure the physical properties of rock and reservoir fluid^{[21,22,23,24]}. Weger et al.^[25] also studied the variation of pore space due to pressure variations by using elastic wave velocity, improving the precision of permeability estimation. By using digital image analysis methods, they produced images of pore space parameters and investigated the petrophysical properties and wave velocities of carbonate rocks. Abdelkarim et al.^[26] analyzed petrophysical properties of carbonate rocks and concluded that the compressive wave velocity was related to pore type and increased with the increase of pressure.

This study makes an attempt to examine the effect of cyclic stress loading on petrophysical properties of carbonate reservoir rock in combination with sonic wave velocities and new approach for CT scan analysis to understand structural alteration in porous media and the storage and transport of reservoir fluids in petroleum reservoirs. To identify heterogeneity, the presence of possible fractures, and changes in porous space, CT scan images and CT number analysis were used as advanced normalized diagrams during loading and unloading processes for the first time. In addition, optical microscopy and scanning electron microscopy (SEM) were also used to analyze changes in crystal structure shape and particle size and degree of compaction in the rock samples during stress loading.

In order to study the effect of pressure on the petrophysical properties and structure of the reservoir rock, two carbonate core plug samples were taken from an oil field in southwest Iran. The selected samples different in lithological composition were taken from two different depths (2974.4 m and 3119.4 m). Since the aim of the experiment was to investigate the effect of pressure loading/unloading on rock samples with different geological structures, the samples were an anhydrite dolomite sample and a calcite sample selected according to macroscopic and microscopic studies. As they were taken from depths adjacent to each other, they have close petrophysical properties. The samples were washed with toluene and methanol using a Soxhlet unit to remove all fluids and salts. Then the plug samples were placed in an oven at 60 °C for 24 hours. A helium porosimeter and air permeameter were used to measure the porosity and permeability of the dried samples, respectively. The results are shown in Table 1.

In order to investigate the changes caused by the loading on the porous media and the structure of the reservoir rock, the sonic wave velocity was measured for the samples. The testing process included several steps that are described in the following sections.

CT scan is a non-destructive radiological imaging technique that measures the attenuation of an x-ray beam as it passes through a material. The incoming radiation may be reflected, scattered, absorbed, and re-emitted as lower energy electromagnetic radiation. It may also be transmitted through the material. The percentage of transmitted radiation is a function of the material thickness, density, and chemical composition^[21,26-33]. The normalized CT number is defined as follows:

where µ is the density of material, dimensionless, µ_w is the density of water, dimensionless, n_CT is the CT number, dimensionless.

Fig. 1 shows the CT scan images of 12 cross sections of the samples before pressure loading. The heterogeneity of and fractures in the selected samples can be observed using the CT scan images. The results of the CT scan indicate that sample A contains a small number of anhydrite nodules, while sample B is relatively homogeneous and free of fractures, in order to better distinguish these nodules, they are outlined in red circles in the figure.

In this study, a CMS-300 apparatus (Core Lab) was used to test the petrophysical properties of the samples. The testing principle of this device is shown in Fig. 2. The device uses a Hassler sleeve to apply radial and axial stress on the sample simultaneously. The range of the applied stress is from 3.5 MPa to 70 MPa. The device uses gas to measure pore volume and permeability^[34]. The measured porosity is effective porosity. The device can automatically correct the air permeability to liquid permeability and the test result is more accurate than that from the conventional steady-state analysis method^[35,36].

In this study, the compression (v_p) and shear (v_s) wave velocities were measured to characterize the elasticity of the reservoir rock. The speed of the compression and shear waves in carbonates are related to the rock matrix, porosity, mineral content, and cementation manner of the grains^[25-26,37]. The geometry of pores, porosity, and sonic wave velocities in carbonate rock, are related to its diagenetic conditions^[38,39]. The sonic wave velocity is an ideal index to measure the compaction of the rock. Since pressure loading affects the factors related with the wave speed, the velocity of these waves can be analyzed to assess the variations resulted from cyclic loading. SonicViewer-SX (Fig. 3) was used to accurately measure the travel time of ultrasonic waves in rock samples to obtain the compression (P) and shear (S) wave velocities at each pressure loading step.

In order to investigate the structure of the sample, the core samples were analyzed by optical microscopy and SEM. Compared with optical microscopy, SEM images have a much larger magnification and higher resolution.

The porosity and permeability changes due to loading of samples A and B are shown in Figs. 4 and 5 and Table 2. It can be seen from Figs. 4 and 5 that the reduction trend in the porosity was the highest in the initial pressure loading stage (step 1 and step 2). Table 2 lists the variations of physical properties of the samples under cyclic pressure loading condition. It can be seen from this table, at the end of the first pressure loading cycle (step 1 and step 2) with pressure loading pressure of 13.79 MPa and 27.58 MPa, the pore volume and permeability of sample B reduced by 5.98% and 31.27% than the original values, respectively, while those of sample A reduced by 4.77% and 4.32%, respectively. In the second pressure loading cycle (step 3 and step 4), the pore volume and permeability of sample B reduced by 4.2% and 13.49% than those at the end of the previous cycle, respectively, while those of sample A decreased by 3.22% and 0.12%, respectively. The pore volume and permeability of sample B continued to decrease as the experiment went on, and the decrease degree of the first two pressure loading cycles was higher than that of the last cycle. For sample A, the pore volume continued to decrease as the experiment continued, but the permeability did not conform to this trend. In the third stage, the permeability increased instead.

The porosity decreased in both samples due to the cyclic loading. When the rock was subjected to loading, it deformed and its compressive strength was reduced. Normally, cyclic loading damages and ultimately fatigues the rock, which directly affects its compressive strength^[14]. During the pressure loading, when passing the yield stress state condition, the sample would take on ductile state, with particles slipping and rotating in the pores. As a result, the grains crush into the weak sorting parts, pores collapse, and the porosity reduces ^[3]. After the initiation of dilatancy under loading, the microcracks also form in the rock sample and increase its permeability^[40,41]. On the other hand, the core sample A had a different lithology and compression of dolomite and anhydrite. Thus, according to their different characteristics, they can show anomalous behavior during loading, which can be found on the boundary between these two lithological compositions and cause a permeable boundary and permeability enhancement.

Wave velocity is a function of the compression and the distance between the rock grains. The denser the rock structure, the higher the wave velocity. The velocity of compression and shear waves of samples A and B under long-term loading are shown in Fig. 6. It can be seen with the increase of loading pressure, micro-fractures gradually close^[42,43,44], as a result, the compression and shear wave velocities increase. Moreover, the acoustic velocity of sample B increased significantly at the first pressure loading step. In the next pressure loading step, the velocity increased at a lower rate. At the initial and final pressure loading steps, the compression wave velocities were 2381.56 and 3360.38 cm/s, respectively, while shear wave velocities were 1380.61 and 2771.31 cm/s, respectively. Since the compression and shear wave velocities increase with the decrease of porosity, when the pore volume of sample B decreased during pressure loading, the porosity variations with the wave velocities are consistent with the results of King’s study results^[45]. At low loading pressures, the wave velocities increased with the increase of hydrostatic pressure due to closer grain contact and closure of pores and microfractures^[46,47]. From the diagram of the sonic wave velocities (Fig. 6), sample A had a variation trend different from sample B. Its compression and shear wave velocities followed irregular trends during the pressure loading. Its compression wave velocity reduced from 3287.5 cm/s to 3114.21 cm/s during the loading stage, while its shear wave velocity increased from 1588.95 cm/s to 1787.61 cm/s.

CT scans were taken from the samples and analyzed as normalized diagrams throughout the sample, the normalization method is shown in formula (1). All cross-section CT images were analyzed, and the minimum, maximum, average, median and mode of the CT number in each image were calculated and plotted (Figs. 7-9). These curves show the statistical analysis results of CT scan numbers for each sample, with the horizontal axis being the normalized length of the sample (0 to 1) and the vertical axis being the CT number (in the Hounsfield unit). As shown in the figures, there is a spectrum of colors from blue to red between the minimum and maximum curves. These colors represent the normalized frequency of the CT numbers from 0 to 1 within the image of each section, which could not be shown in the 2D plot unless a color scheme is used. At low frequencies (usually near the minimum and maximum values), the color of the plot is blue, whereas at the high frequencies (as the mode reaches a maximum and equals to 1), the color of the plot is red. The parameters mentioned above and the color spectrum show the heterogeneity and fracture patterns of the samples.

Figs. 7 and 8 show the normalized CT scan histogram of samples A and B at all the pressure loading steps. In the figures, the red lines represent the maximum CT number and the blue lines are the minimum CT number in the cross-sections (the maximum CT number of the CT scan device was set at 3071 HU). The black line represents the average CT number along the sample axis. The mode section shown in green represents the frequency of the CT number per slice in the measured sample. For example, as shown in the graphs, in some cross-sections the green section is below the average value, which indicates a lower frequency of the CT number in the desired slice. Also, there are different spectra in the center between 0 (blue) and 1 (red). It should be noted that overlapping of the mode, median, and average on a line represents the homogeneity of the sample.

As shown in Fig. 7, the variation of the CT number of sample A indicates significant changes during the various pressure loading steps. Prior to the third loading step, there was no significant change in the CT number, but during the third loading step, an intense peak was observed at the minimum and maximum value of the CT number, indicating that fracturing occurred inside the sample. This diagram also confirms that the structure of sample A changed abnormally and the permeability increased suddenly after the second step of cyclic pressure loading.

SEM was used to study the quality of the pore structure, the existing porosity, as well as investigate the microfractures and their expansion in the samples. Figs. 10 and 11 show the SEM images of core samples before and after pressure loading. In Fig. 10, the changes in structure and porosity of samples A and B can be observed. Before pressure loading, Sample A had dissolution pores (indicated by yellow lines) and angular pores (indicated by blue lines, see Fig. 10a). After loading pressure, the pore space decreased and porosity reduced (Fig. 10b-10e). In the calcite sample B, phenomena of pore volume reduction, particle crush, and changes of pore size and grains can be observed during the pressure loading process (Fig. 11).

As shown in Fig. 4 for sample A, the permeability showed no significant decrease despite the relatively significant porosity reduction at the initial pressure step at 13.79 MPa. The permeability trend of sample A in this relatively low-pressure step was a function of porosity. This indicates that for sample A at the initial step at 13.79 MPa, the porosity reduced due to pore collapse and deformation of the pore spaces. The permeability did not reduce much because the grain crushing did not cause blockage of the microcracks. On the other hand, it can be claimed that by increasing the pressure to 27.58 MPa, the closure stress of the fractures could not be overcome.

It should be noted that during the transition from step 2 to step 3 in sample A, the permeability increased even though the porosity reduced for all steps. This indicates that despite collapsed pores and a possible slight decrease in the permeability in some places, the overall permeability increased due to the creation of major fractures or the joining of connected microcracks that created a new conduit. This phenomenon can be observed in SEM images of sample A in Fig. 12. Based on the changes in permeability (Table 2), during the third step of pressure loading (Fig. 13d), the increase in permeability was due to the opening of fractures and the coupling of paths connecting pores. This conclusion is consistent with the change of permeability. These changes have also been interpreted by processing the microscope images of the dolomite-anhydrite boundary shown in Fig. 13. By comparing the changes of SEM images before pressure loading and after various pres-sure loading steps on the boundary between dolomite and anhydrite, it can be seen that some changes occurred. With the increase of loading pressure, fractures formed and expanded at the boundary, so came the permeable boundaries. Consequently, during the third stage of loading, the crack opening and permeability increased. Also, the reduction in the pore volume by loading and applying pressure is shown in Fig. 10, and was due to the closing of vuggy (yellow arrows) and angular (blue arrows) porosities.

At step 5, the permeability of sample A was approximately equal to that in step 3, indicating that the fractures formed in step 4 were blocked. As can be seen in the SEM image in Fig. 12(e), the decrease in the permeability was a result of the closure of fractures and the placement of crushed particles in the crack path. During step 6, despite the unloading, the permeability did not increase and was almost constant, which was due to the absolute plugging of these new fractures induced by the collapse and grain slippage in the pore paths. In sample A, from step 2 to step 3, although the fractures closed, the permeability restored to its initial value when unloading. At step 5, the fractures were closed with the increase of pressure, although the porosity from step 5 to 6 did not reduce noticeably, and the permeability almost remained constant. The cyclic pressure loading and unloading caused shearing and collapse in the grains and irreversibility of the permeability and minor alterations in porosity. This crushing resulted in the disruption of the fractures that re-open the pathway for permeability during unloading.

It can be seen from Fig. 5, in sample B from step 2 to 3, the porosity decreased and the permeability was almost constant during the pressure unloading. This proves that the collapsed pore spaces were not in the mainstream of connected pores, which is why the conductivity of the system did not change. However, maybe new microcracks were generated in the system, which increased the permeability and compensated for the permeability reduction resulted from the pore collapse. The crushing of pores and pore size changes can be observed in the SEM images in Fig. 11a to 11e.

The variation in the velocity of the sonic waves, in particular the compressive velocity of sample A, was due to the changes in the pore geometry and pore spaces where vuggy sections were present (Fig. 10). Also, due to the lithological and mineral composition of sample A and the difference in the density of dolomite and anhydrite and their compression properties, different values ​​of velocities under different loading conditions were observed. Applying pressure to the sample created two mechanisms that affected the velocity of the elastic waves. These two mechanisms are the generation of microcracks and closure of available cracks. For example, in sample A, an extraordinary increase in the sound velocity was observed from step 4 to step 5, which was due to the complete closure of the fractures. In sample B, the V_p did not change significantly from steps 5 to 6, but V_s increased dramatically.

In sample B, during the loading step 1 to 2, the closure stress of the microfractures was smaller than for sample A. Upon increasing the pressure to 13.79 MPa, considering the two mechanisms of deformation of the porous media as well as the closure of the microcracks, the velocity of the waves dramatically increased. In the loading step 2 to 3, there was a sharp decrease in porosity and permeability and no change in the wave velocities. In step 1 to 2, only the elastic deformation of the pores and closure of the microcracks occurred. However, after increasing the pressure to 27.58 MPa, the elastic stage passed and caused pore collapse and grain crushing despite the porosity and permeability reduction. This inelastic deformation has a destructive effect on the velocity of elastic waves, such as P and S^[48]. In this case, due to grain crushing, pore collapse, and grain slipping, the porosity was reduced and consequently the wave velocities increased. Also, due to grain crushing, the growth of microcracks, and the increase in fracture density, the wave velocities decreased.

At the initial pressure step up to 13.79 MPa for sample A, the fracture density did not reduce, which consequently did not affect the permeability, unlike for sample B. However, the sound wave velocity decreased to a noticeable amount because crushing in the pores occurred. Upon increasing the pressure up to 27.58 MPa in sample A, the sound wave velocity was almost constant due to the reduction of the fracture density mechanism and grain crushing. Upon decreasing the load to 13.79 MPa, the wave velocity decreased due to the opening of the cracks and the growth of microcracks.

Fig. 9 shows the average value of the CT number for samples A and B based on the normalized length of the plugs at different pressure steps. The difference in mean CT numbers in sample A is more in the first half of the sample (along 3 cm), especially during the second loading stage, which confirms the increase in the permeability due to the expansion of the microcracks. Also, most changes that occurred in the first half of the sample length were due to the existence of anhydrite in this section, as shown in Fig. 1 and SEM images in Fig. 13. This means there was a higher heterogeneity in sample A and more anomalous alteration in the permeability due to loading. For sample B, upon loading in the first stage (13.79 MPa), the average CT number increased due to the closure of the pores and the compression of grains. In the second stage (27.58 MPa), the CT number decreased in most sections, by creating microcracks. In the next step, incremental growth was observed. In general, the proximity of changes in the various loading steps represents the homogeneity and the higher compressibility of the sample.

The changes in petrophysical and structural properties of two carbonate reservoir samples with different lithological compositions under hydrostatic cyclic loading were described as follows: (1) The identical behavior of the samples in reducing pore volume due to cyclic loading was observed, which can be due to the slippage and rotation of the particles in the pores, as well as the crushing of the particles and pore collapse. (2) Changing and decreasing permeability was observed. Thus, for sample A, which was composed of dolomite with anhydrite nodules, this was observed in the third loading stage due to the start of compression. This was especially case because of the presence of two lithological combinations in the rock sample with different compression properties. As a result, permeable boundaries were formed at the boundary between two compounds and cracks were induced that ultimately increased the permeability. (3) The changes in the compression and shear wave velocities under cyclic loading indicated that the variation in pore geometry, porosity type, and especially the heterogeneity in the lithology composition of the reservoir, can be observed and evaluated using these waves. Moreover, the effect of loading on elastic waves was evident based on the formation of microfractures, their closure, and compression of the sample. (4) Analysis of CT scans and diagrams showed heterogeneity and the presence of fractures during various loading steps. Hence, when increasing the permeability during loading at the third stage, the CT scans showed severe peaks at the maximum and minimum CT numbers with a sharp decrease in the CT. Therefore, the CT is a good tool for investigating the structural changes in porous media. (5) The analysis of microscopic and SEM images revealed the structural changes of the tested samples, including porosity variations in the rock matrix crushing due to cyclic hydrostatic loading. In particular, the growth of the microfractures during loading at the third stage and the closure of fractures when loading the fourth stage were observed. Similarly, the microscopic investigation of the presence of a permeable boundary on the dolomite and anhydrite boundary proved the hypothesis that samples with different lithological compositions cause fluctuations in the petrophysical properties. (6) There was a good convergence for all the measurements, including porosity, permeability, sonic wave velocities, SEM images, CT scan, and optical microscope images.