Variable secondary porosity modeling of carbonate rocks based on μ-CT images

2.1 Dilation and erosion in mathematical morphology

Mathematical morphology is a standard method used in image processing. The principle of mathematical morphology method is to use a structure element with a specific shape to measure and pick up the equal shape in the image to analyze the image. There are four basic operation modes: dilation, erosion, opening, and closing [26]. The opening and closing operations have been used on sandstones to generate fluid distributions in the pore space [27].

In mathematical morphology, dilation operation and erosion operation are opposite processes. The dilation operation enlarges the target image, while the erosion operation reduces the target image. The mathematical morphology is based on set theory. Sets in mathematical morphology represent the shapes which are manifested on binary or gray tone images. A and B are two sets; A ⊕ B means the dilation of A by B, while A ⊖ B means the erosion of A by B. Eq. 1 and 2 are expressions of dilation and erosion,

respectively.

In Eq. 1 and 2, B_x is the translations of set B by x; ∅ is null. Figure 1 is a 2D sketch of dilation and erosion. Figure 1(a) shows an original image M. In this case, B is a disc with a fixed radius, and the origin is the center of the circle. After the dilation of M by B, the image M⊕ B is as the grey part shown in Figure 1(b) The original image is enlarged with a certain radius outward. Figure 1(c) is another original image N. The grey part shown in Figure 1(d) is N ⊖ B, i.e the erosion of image N by B. In this case, M and N are images, and B is structuring element. In our study, we use a sphere structuring element for 3D processing.

Figure 1

Sketch of dilation and erosion

2.2 Fractional Brownian motion

Experimental evidence [28] shows that the roughness of natural fracture internal surfaces in rocks follows self-affine fractal statistics. In this paper, we adopt the idea of considering a rock surface without overhangs. The surface height is described by a function h(x, y), with the coordinates x and y lying in the mean plane of the fracture [10]. Self-affinity of the rough surfaces implies scale invariance. When x → λ_xx, y → λ_yy, h → λ_hh, we assume that λ_x = λ_y = λ and λ_h = λ_H, and then we have the following equation:

where λ is the scale factor, H is the roughness (Hurst exponent). For the surfaces of fractures presented in rocks, the roughness exponent is around H = 0.8, which has no relationship with the lithology and fracture mode [29].

Fractional Brownian motion (FBM) is a standard mathematical model for self-affine fractals, which fits the following equation:

where G_H (X) is a random walking which satisfies the Fractional Brownian motion; E is mathematical expectation; His the roughness exponent; σ is standard deviation; σ⁰ is a random initial standard deviation between 0~1; j is the times of iterative operations.

According to central-limit theorem, Eq. (4) can be written as:

standard deviation σ can be calculated as follow equation:

There are many algorithms to generate FBM data. In this paper, we adopt a corrected successive random addition (CSRA) algorithm to generate the FBM surfaces [30, 31].

The workflow of this algorithm for 2D square surface is as follows (Figure 2):

Figure 2

Point denoting and value assignment order of SRA

1. First select random numbers from and add them to the 4 corner points of the square, denoted by number 1, where stands for a Gaussian random number generator with a mean of 0 and a variance σ02.

2. Calculate the midpoint denoted by number 2 by linear interpolation based on the 4 corner point values.

3. Select random numbers from (σ12can be derived from Eq. (6)) and add them to all the existing points;

4. Calculate the midpoints on each side of the square denoted by number 3 by linear interpolation using two adjacent corner values obtained in step 3;

5. Select random numbers from N(0σ22)(σ22can be derived from Eq. (6)) and add them to all existing points.

6. Repeat the process, keep performing linear interpolations and adding random numbers from N(0σj2)up to the nth level, resulting in a total of (2ⁿ+1)² points. Then the generation of the 2D rough surface is done (e.g, Figure 3).

3.1 μ-CT based 3D digital rock

In this study, we use a 3D image of a carbonate rock (named C1) obtained using Microcomputed X-ray tomography (μ-CT) from Imperial College London [32]. The segmented 3D carbonate digital image is shown in Figure 4, which has 400×400×400 voxels. The resolution is 2.85μm/voxel, and the porosity is 17.6%.

Figure 3

Rough surface created by using SRA algorithm

Figure 4

The μ-CT scanned digital rock sample C1 of carbonate rock used in this study

3.2 Modeling of carbonate digital rock with different porosities

One most important component of carbonate reservoir pore space is secondary porosity formed by corrosion. Corrosion process in the carbonate rocks is similar to the erosion operation in mathematical morphology. Therefore, we use erosion to mimic the later stages in a corrosion process. Based on the digital rock model obtained by μ-CT scanning and using the erosion operation, we generated the digital carbonate rocks with different porosities.

Figure 5 shows an example result of digital carbonate rocks with different porosities after the erosion. In these digital rocks, pixels with value 0 (shown as black parts)

Figure 5

Digital rocks and their slices of carbonates with different porosities generated after the erosion process. (a) Original image (Φ= 23.26%); (b) Image after 1 erosion (Φ = 29.22%); (c) Image after 2 erosion (Φ= 34.93%); (d) Image after 3 erosion (Φ = 40.4%)

represent pore space, while those with value 1 (shown as lighter parts) represent the matrix. The structuring element is a sphere with a radius of 1 pixel. Figure 5a is the original core and its one slice, and Figure 5b~d is the cores and their slices after 1~3 times of erosion. As we can see in Figure 5, after the sphere goes around the pore space, the matrix is eroded, resulting in an increase of porosity. And during the process, pores used to be close to each other can merge as one larger pore. Through this way, we can increase the porosity of the carbonate digital rock models.

With the erosion, we generated digital rocks with higher porosities. Meanwhile, by using the dilation process to dilate the matrix, digital rocks with lower porosities can be acquired. The structuring element is the same as in the erosion process. Figure 6 shows the dilation results. Figure 6a is the original core and its one slice, and Figure 6b~d is the cores and their slices after 1~3 times of dilation. As we can see in Figure 6, after the sphere goes around the pore space, the matrix is dilated, resulting in a decrease of porosity. Some small pore space vanishes into the matrix. This process can be seen as reverse of a corrosion process. The cores with lower porosities can be seen as earlier stages of the corrosion process. Through thisway,

Figure 6

Digital rocks and their slices of carbonates with different porosities generated afterdilation process. (a) Original image (Φ= 23.26%); (b) Image after 1 dilation (Φ= 17.76%); (c) Image after 2 dilation (Φ= 13.5%); (d) Image after 3 dilation (Φ= 10.22%)

we can decrease the porosity of the carbonate digital rock models.

3.3 Modeling of fractured rocks

Fractures are another critical space for carbonate reservoirs. They can be both reservoir space and seepage channels. But due to the difficulty in coring and doing experiments for fractured formations, the research about fractured carbonate reservoirs hasn’t been well done before. Digital rock technology can be introduced to provide fractured rock models whichcan be the basis of numerical simulations. In this study, we use fBm to generate fractures with different widths and angles.

For a digital rock with 400×400×400 voxels, we choose a 400×400 pixels area in the rough surface z=f (x,y) generated by using FBM in Section 2. Then we set it as the lower surface of the fracture. Set the width of the fracture as d, then the upper surface Z=z+d. Because 0 represents pore space, we set the pixels between upper and lower surfaces as 0. By changing the value of d, we can get different fracture widths. By setting the pixels where the fractures are supposed to be as 0 in the original image, we can obtain the digital rock models with fractures different in widths (e.g, Figure 7).

Figure 7

Sections of fractured carbonate digital rock models with different fracture widths

According to the attitudes of the fractures, they can be divided into “high angle” (dip angle > 75^∘), “bevel” (15^∘ < dip angle < 75^∘) and “low angle” (dip angle < 15^∘) fractures. To mimic different attitudes, by rotating the fractures, we can obtain fractured digital rocks with varying angles of dip (e.g, Figure 8).

Figure 8

Models of fractured carbonate digital rock models with different fracture angles

4.1 Pore geometry, numbers and porosity changes of models after erosion and dilation

Because the secondary corrosion pore space in carbonate rocks is formed through the corrosion process which is similar to dilation and erosion operations in mathematical morphology, the dilation and erosion processes are adopted in this research to simulate different stages of eroded pores. Our goal is to have variable porosity but retain a similar pore structure. Therefore, we use the maximal ball algorithm [32] to extract the pore network (Figure 9) and statistic the pore and throat numbers of all radiuses. The balls and sticks in Figure 9 represent pores and throats, respectively, and different colors show different pore sizes. From blue to red, the pore size increases.

Figure 9

Pore network of C1 digital rock model

From the network models, the pore/throat radius distribution curves of the digital rock models with different porosities can be obtained (Figure 10). In Figure 10, the frequency in y-axis stands for the count of the radius in x-axis divided by the whole count. C11~C15 are the names of digital rock models after erosion, and C1_1~C1_7 are the ones after dilation. As we can see in Figure 10, the peak of the curve slightly shifts right in the erosion process, while it moves left in the dilation process. This reveals the pores are somewhat enlarged or shrunken during the respective operations.

Figure 10

Pore/throat radius distribution curves of the digital rock models

The numbers of pores are counted (Figure 11). In Figure 11, each dot represents a dilation or erosion processing point. In both processes of erosion and dilation, the pore numbers decrease. This also indicates that in the dilation process, the small pores vanish; and in the erosion process, the nearby pores merge into bigger ones. Figure 12 shows the porosity change rate during the operations. As we can see in Figure 12, the more the processes have been done, the slower the erosion and dilation process is. Because each process only adds on or peels off a fixed thickness of one voxel on the pore surface, the speed of change relates only to the surface area. As a result of small pores vanishing or nearby pores merging, the pore surface decreases in both processes.

Figure 11

Pore numbers of digital rocks with different porosities

Figure 12

Porosity change variation with the processing times

4.2 Porosity-permeability relation

The relationship between porosity and permeability is always concerned by petroleum engineers and petrophysicists. Many studies reveal that in on a particular area or for one specific kind of carbonate, the permeability can be a power function of porosity (e.g, Rashid et al., 2015). We have some carbonate rock samples with laboratory experimental permeability results (Figure 13). The trend is also power exponential. As we can see, there are rock samples with minimal porosity and high permeability. This reveals there might be micro-fractures in those samples. Based on the network models, the permeabilities of the digital rock models are calculated by using the single-phase flow permeability calculation method based on Darcy’s law [34]. The simulation results are shown in Figure 14. As we can see in Figure 14, the trend is similar to the experimental results of carbonate rock samples from a specific area acquired in the laboratory. Because the simulations are based on the same sample, the trend of the simulation results is more evident than the laboratory results. This indicates the method used to modify the porosity in this paper is reasonable. The models can be used in the next-step physical simulations.

Figure 13

Laboratory experimental permeability results of carbonate rock samples

Figure 14

Permeability simulation results of the models generated after erosion and dilation

4.3 The effect of B’s radius on porosity changing

In this study, we adopt a sphere with a radius of 1 voxel as the structuring element B. The porosity changes by one layer of voxel on the pore surface during each process. Because the resolution of this 3D digital rock is 2.85μm/voxel, the pores’ radiuses change 2.85μm(by add or minus) each time. In a pixel image, processing with the sphere structuring elements with different radius can cause different results. We compare the dilation and erosion results of the digital rock by a sphere with a radius of one voxel for five times and by a sphere with a radius of five voxels for one time. The statistic results are shown in Table 1. C13 and C1_3 are the results of erosion and dilation by B with a

radius of one voxel for three times, and C13r and C1_3r are the ones by B with a radius of three voxels for one time. There are differences in porosity, permeability, and pore numbers. The porosity and permeability change more

when B’s radius is three voxels than one voxel. Figure 15 shows the pore networks of the four digital rocks. From the network images, we can see some differences in pore connection. In the dilation process, when B’s radius is three voxels, the smaller throats and pores are less. This makes

Figure 15

Pore networks of digital rock models with different porosities

the permeability lower. This indicates when a large radius is adopted for a sphere structuring element, more details will be erased earlier than using a one-voxel-radius structuring element.

4.4 Method advantage and limitation

Based on the μ-CT image, erosion and dilation operations are adopted to generate different porosities in this study. Other mathematical methods have been used to construct models with different porosities before, including stochastic methods [18, 19], process-based methods [20]. When using these methods, each time of the process, the porosities are different, and so are many other parameters. The pore shape of one model can’t be guaranteed to maintain the same way as the others. In this study, operations are based on the same μ-CT image. Therefore, the basic pore structures and distributions of the models are the same. And based on the dissolution process, many researchers have done the numerical study of carbonate dissolution by CO2 and its effect on the transport properties [21, 22, 23, 24, 25]. These process-based methods can provide us accurate models with different porosities after the present dissolution stage. In this study, the erosion operation can approximately simulate the corrosion effect, while the dilation operations can restore the matrix before the current status. By using dilation and erosion, both higher and lower porosities can be generated.

However, dilation and erosion operations are used under assumptions, namely: (a) the rock corrosion happens evenly and simultaneously at each direction on each pore surface, and (b) the pore will be corroded no matter it is connected to the primary fluid channel or not. Therefore, these methods have nothing to do with physical properties such as the heterogeneous of minerals and gravity. More studies need to be done in the future.

4.5 Reliability of fractured models

We use the FBM method to model the fractures in carbonate rocks based on μ-CT images. The surfaces of natural fractures in rocks have very high self-correlation and self-affine property. The FBM is a standard self-affine model. By using this method, we can generate a fracture-like rough surface. Zhao et al. (2013) [10] applied the FBM method to create fractures in sandstones. And based on the fractured models, the electrical simulations were done. By comparing the simulation results with the experimental results, the FBM method was proved efficient. In this study, we adopt this same method into the carbonate rocks. Because the method is based on mathematics, and it has nothing to do with the rock properties, the usage of the FBM method in carbonate rocks is also workable. The fractured digital rock models can be used in the next step physical simulations.

μ-CT is one common but not the only way to build 3D digital rocks. Because the methods and workflows introduced in this paper are purely based on mathematics, they can also be used in carbonate digital rocks reconstructed in other ways, such as using process-based methods.