Pore-Scale Characterization and PNM Simulations of Multiphase Flow in Carbonate Rocks

: This paper deals with pore-scale two-phase ﬂow simulations in carbonate rock using the pore network method (PNM). This method was used to determine the rock and ﬂow properties of three different rock samples, such as porosity, capillary pressure, absolute permeabilities, and oil–water relative permeabilities. The pore network method was further used to determine the properties of rock matrices, such as pore size distribution, topological structure, aspect ratio, pore throat shape factor, connected porosity, total porosity, and absolute permeability. The predicted simulation for the network-connected porosity, total porosity, and absolute permeability agree well with those measured experimentally when the image resolution is appropriate to resolve the relevant pore and throat sizes. This paper also explores the effect of the wettability and fraction of oil-wet pores on relative permeabilities, both in uniform and mixed wet systems.


Introduction
In the context of reservoir characterization, identifying the pore structure and understanding flow and transport processes in sub-surface geological rocks are nontrivial tasks. For instance, relative permeability, which plays a key role in reservoir characterization and production, is measured in the laboratory using either steady-or unsteady-state methods, which are lengthy and costly. These experimental methods involve injecting different phases into core samples and measuring flow rates [1][2][3][4][5]. Measuring relative permeability in the laboratory also requires a large number of high-quality core samples from a targeted reservoir [5,6].
Nowadays Digital Rock Physics (DRP), applied to reservoir rock characterization, is becoming a prominent method to determine physical and flow properties [7][8][9][10]. DRP relies on high-resolution three-dimensional digital images of pores and grains structures, obtained using X-ray micro-computed tomography (micro-CT). These structures are used to conduct multiphase flow numerical simulation at the pore scale and to estimate rock properties such as relative and absolute permeabilities. The success of the DRP approach is related to the development of advanced imaging techniques and computational capability. Micro-CT, which is a non-destructive imaging technique, can resolve pore structures of a few microns [5,11].
The numerical simulations of the multiphase flow through pore structures are of two types, namely, direct simulations and pore network methods. Direct simulation methods include the lattice Boltzmann method (LBM) and the finite volume method (FVM). These methods solve the multiphase flow within pore structures but require considerable computational resources and are time-consuming. These drawbacks grow with the resolu-tion of a pore structure. Pore-scale direct simulations of multiphase flows are costly and time-consuming [11][12][13].
The pore network method (PNM) is a computationally cost-effective method to predict multiphase flow and transport properties. The pore network method was introduced first by Fatt [14], based on regular tubes and capillary pressure. With time, the pore network method became more sophisticated to accommodate irregular and complex pore networks in real rock. Different network extraction algorithms such as the medial axis algorithm [15,16], the maximal ball algorithm [17,18], and the SNOW algorithm [19] have been developed to generate the equivalent structure of porous media. The differences between these algorithms are based on different pore space characterization and pore size assignment methods. The output pore network is then used to conduct flow simulations and to investigate different processes such as wettability, multiphase flow, reactive transport, and no-Darcy flow [20][21][22][23][24][25].
One can find in the literature many investigations that have used PNM to conduct various multiphase flow simulations in porous media, to predict the permeability, and to investigate the effect of various factors such as wettability and pore connectivities on the permeability. For instance, Valvatne and Blunt used a network representation of processbased simulation Berea sandstone to predict the relative permeability for different samples with different wettabilities [26]. PNM was used to predict two-phase water-wet and threephase gas injection relative permeabilities in Berea sandstone samples [27]. Al-Kharusi and Blunt used a PNM on a statistically representative 3D image from a 2D thin-section carbonate image to predict the relative permeabilities of heterogeneous carbonate rocks [28]. Zahaf et al. predicted the relative permeability with networks extracted from 3D images of real carbonate samples [29]. A pore network obtained from real 3D carbonate images was used by Gharbi and Blunt to investigate the effects of wettability and connectivity on relative permeability in mixed wet systems [30]. Roth et al. used pore networks extracted from 54 samples and studied the effects of rock types and the porosity on water-oil relative permeability in the primary drainage process [31].
There are few papers where PNM has been used to predict permeability in both uniform and mixed wet systems. Moreover, to the best of the authors' knowledge, most of these investigations used PNM obtained from low-resolution images. For instance, Zhao et al. and Valvatne and Blunt used PNM to estimate the relative permeability in both uniform and mixed wet systems [26,32]. In Zhao et al.'s study [32], the resolution of the images ranged from 5.345 to 9.996 µm; with such low resolution, nano-and micropores and throats, which are prominent in carbonates, cannot be completely resolved. This means the networks extracted in [32] cannot be representative of their carbonate rock. The resolution of the throat is crucial to accurately predict capillary pressure. Failing to resolve nanoand micropores and throats should result in an underestimation of the capillary pressure curve. For instance, Valvatne and Blunt had to adjust the pore size distribution to match the capillary pressure [26]. Piri and Blunt adjusted the contact angle to match the experimental relative permeability to obtain good agreement between the results of the numerical simulation and the experiment [27]. Zahaf et al. matched the measured capillary pressure by adjusting the contact angle [29]. In the present study, we used a much higher image resolution (0.94 µm and 0.81 µm); this resolved micropore and microthroat in the carbonate rocks. Hence, we expect our extracted pore networks to be more representative of the rock samples in hand. Moreover, the capillary pressure measured in the laboratory should be predicted without adjusting any parameter such as pore size or contact angle. Table 1 highlights the novelty in our present study.
In this study, we conducted pore-scale two-phase flow simulations using PNM models of carbonate rocks. The carbonate core samples from oil and gas reservoirs presented different heterogeneities. First, we used pore network models extracted from 3D micro-CT images to analyze rock structure properties such as pore size distribution and topological structure. Then, we conducted simulations of multiphase flows through the pore network models. We were able to predict the porosity, absolute permeability, and capillary Energies 2021, 14, 6897 3 of 20 pressure curve. These results were validated using experimental data. We also predicted the oil-water relative permeability and analyzed the effect of wettability on the relative permeability in both the uniform and mixed wet systems.

Rock Samples and Image Acquisition
In this study, we selected three different carbonate rock samples, including Silurian Dolomite (SD), Albion-4 carbonate rock (A-4), and the TOTAL carbonate rock (TC), all from TOTAL reservoirs (see Figure 1). The selected rock samples were trimmed to fit a diameter of 1.5 inches. Smaller plugs were then extracted from the individual samples and scanned using Zeiss Xradia Versa 500 Micro-CT, Oberkochen, Germany, at 5.32 µm for SD, 0.81 µm for A-4, and 0.94 µm for TC, respectively (see Figure 2). CT images to analyze rock structure properties such as pore size distribution and topological structure. Then, we conducted simulations of multiphase flows through the pore network models. We were able to predict the porosity, absolute permeability, and capillary pressure curve. These results were validated using experimental data. We also predicted the oil-water relative permeability and analyzed the effect of wettability on the relative permeability in both the uniform and mixed wet systems.

Rock Samples and Image Acquisition
In this study, we selected three different carbonate rock samples, including Silurian Dolomite (SD), Albion-4 carbonate rock (A-4), and the TOTAL carbonate rock (TC), all from TOTAL reservoirs (see Figure 1). The selected rock samples were trimmed to fit a diameter of 1.5 inches. Smaller plugs were then extracted from the individual samples and scanned using Zeiss Xradia Versa 500 Micro-CT, Oberkochen, Germany, at 5.32 μm for SD, 0.81 μm for A-4, and 0.94 μm for TC, respectively (see Figure 2).

Laboratory Measurements
Laboratory measurements of porosity and absolute permeability were obtained for each of the samples as follows. Rock porosity was determined using a helium porosimeter based on Boyle's law, while the permeability was measured by injecting air. The trimmed rock samples were then subjected to mercury injection capillary pressure (MICP) tests to evaluate the connected porosity, pore size distribution, normalized permeability, and capillary pressure as a function of the mercury saturation, as shown in Figure 3. Table 2 summarizes the experimental results obtained for each carbonate sample.
The pore throat size distribution plays a key role in determining the rock's normalized permeability curves, as shown in Figure 3a,c,e, which were determined from the pore throat radius similar to the capillary pressure curves. The normalized permeability curves indicate that only pore throats of 1 μm in size and above contribute to the porous medium permeability. However, Figure 3a shows that the pore throat size distribution for the SD sample ranged from 1 μm to approximately 35 μm, with a prevalence of pores throats of a 10 μm radius, while Figure 3c indicates that the A-4 sample has relatively small pores throats distributed from 0.01 to 5 μm, with several peaks centered around 0.8 μm. Figure  3e features the TC sample pore throat size distribution ranging from 0.001 to 40 μm, which includes more nanopores, more peaks, and a wider distribution. These features are very important in determining the flow and elastic properties of the rock samples. The pore throat size distribution and the normalized permeability curves indicate what image resolution is required in this study. A resolution of 5.32 μm for the SD sample can resolve most of the pore spaces. However, a sub-micron resolution is required for the A-4 and TC samples. Furthermore, in our previous digital rock physics studies, we showed that the nanopores do not contribute to the flow properties of the rock samples [7][8][9]. Therefore, resolving nanopores or including them in the simulation will not affect the estimations of the flow properties such as absolute and relative permeability. However, resolving nanopores and including them in the simulation is important for the prediction of the rock's elastic properties [10,33]. It was shown that the inclusion of nanopores and whether they are empty, filled with air, water, or oil afford different rock elastic properties. This is very important for the interpretation of seismic data in reservoir exploration and characterization. Therefore, the used micro-CT imaging resolutions of 5.32 μm for SD, 0.81 μm for A-4, and 0.94 μm for TC samples are sufficient to extract accurate pore networks of the three samples for our study of multiphase flow properties. As highlighted in Table 1, such high resolutions have not been used in previous works found in the literature.

Laboratory Measurements
Laboratory measurements of porosity and absolute permeability were obtained for each of the samples as follows. Rock porosity was determined using a helium porosimeter based on Boyle's law, while the permeability was measured by injecting air. The trimmed rock samples were then subjected to mercury injection capillary pressure (MICP) tests to evaluate the connected porosity, pore size distribution, normalized permeability, and capillary pressure as a function of the mercury saturation, as shown in Figure 3. Table 2 summarizes the experimental results obtained for each carbonate sample.
The pore throat size distribution plays a key role in determining the rock's normalized permeability curves, as shown in Figure 3a,c,e, which were determined from the pore throat radius similar to the capillary pressure curves. The normalized permeability curves indicate that only pore throats of 1 µm in size and above contribute to the porous medium permeability. However, Figure 3a shows that the pore throat size distribution for the SD sample ranged from 1 µm to approximately 35 µm, with a prevalence of pores throats of a 10 µm radius, while Figure 3c indicates that the A-4 sample has relatively small pores throats distributed from 0.01 to 5 µm, with several peaks centered around 0.8 µm. Figure 3e features the TC sample pore throat size distribution ranging from 0.001 to 40 µm, which includes more nanopores, more peaks, and a wider distribution. These features are very important in determining the flow and elastic properties of the rock samples. The pore throat size distribution and the normalized permeability curves indicate what image resolution is required in this study. A resolution of 5.32 µm for the SD sample can resolve most of the pore spaces. However, a sub-micron resolution is required for the A-4 and TC samples. Furthermore, in our previous digital rock physics studies, we showed that the nanopores do not contribute to the flow properties of the rock samples [7][8][9]. Therefore, resolving nanopores or including them in the simulation will not affect the estimations of the flow properties such as absolute and relative permeability. However, resolving nanopores and including them in the simulation is important for the prediction of the rock's elastic properties [10,33]. It was shown that the inclusion of nanopores and whether they are empty, filled with air, water, or oil afford different rock elastic properties. This is very important for the interpretation of seismic data in reservoir exploration and characterization. Therefore, the used micro-CT imaging resolutions of 5.32 µm for SD, 0.81 µm for A-4, and 0.94 µm for TC samples are sufficient to extract accurate pore networks of the three samples for our study of multiphase flow properties. As highlighted in Table 1, such high resolutions have not been used in previous works found in the literature.

Image Processing
The image processing protocol includes image denoising, removal of artifacts, and classifying pixels into representative clusters [33]. Several image segmentation techniques can be found in the literature. These techniques can be classified into either manual or automatic segmentation methods. The manual segmentation algorithms depend on the operator's observation and experience. The results of the manual segmentation may be optimal for one sample and may not be for another [34]. On the contrary, automatic segmentation algorithms are more efficient and less subjective. In this study, we used the Interactive Overly Threshold module in the PerGeos software package to segment the images into pores and solid grains. Figure 4a-c represent 2D raw images of 500 × 500 pixels of the SD, A-4, and TC samples. Figure 4d-f represent the associated segmented images. The 2D segmented images cannot reveal connected and unconnected pores. Therefore, we applied the same segmentation algorithm (Interactive Overly Threshold module in Per-Geos) to 3D images of 500 3 voxel cubes extracted from full-size images of the three samples. Figure 5 represents the 3D volume rendering of the samples. The gray color in the first column represents the rock matrix, while the blue color represents the total pore volume. In the second column, the gray and the blue colors represent the unconnected and connected pores, respectively. The connected pores, which represent the flow channels

Image Processing
The image processing protocol includes image denoising, removal of artifacts, and classifying pixels into representative clusters [33]. Several image segmentation techniques can be found in the literature. These techniques can be classified into either manual or automatic segmentation methods. The manual segmentation algorithms depend on the operator's observation and experience. The results of the manual segmentation may be optimal for one sample and may not be for another [34]. On the contrary, automatic segmentation algorithms are more efficient and less subjective. In this study, we used the Interactive Overly Threshold module in the PerGeos software package to segment the images into pores and solid grains. Figure 4a-c represent 2D raw images of 500 × 500 pixels of the SD, A-4, and TC samples. Figure 4d-f represent the associated segmented images. The 2D segmented images cannot reveal connected and unconnected pores. Therefore, we applied the same segmentation algorithm (Interactive Overly Threshold module in PerGeos) to 3D images of 500 3 voxel cubes extracted from full-size images of the three samples. Figure 5 represents the 3D volume rendering of the samples. The gray color in the first column represents the rock matrix, while the blue color represents the total pore volume. In the second column, the gray and the blue colors represent the unconnected and connected pores, respectively. The connected pores, which represent the flow channels relevant to the present multiphase flow study, are shown in the third column. This last column shows that the flow channel volume increases from SD to A-4 to TC, respectively. This is further shown by the connected porosity values in Table 3, which are approximately 14%, 19%, and 25%, respectively.
Energies 2021, 14, x FOR PEER REVIEW 6 of 20 relevant to the present multiphase flow study, are shown in the third column. This last column shows that the flow channel volume increases from SD to A-4 to TC, respectively. This is further shown by the connected porosity values in Table 3, which are approximately 14%, 19%, and 25%, respectively.  relevant to the present multiphase flow study, are shown in the third column. This last column shows that the flow channel volume increases from SD to A-4 to TC, respectively. This is further shown by the connected porosity values in Table 3, which are approximately 14%, 19%, and 25%, respectively.

Pore Network Model
We used the pore-based algorithm in the PerGeos software to extract pore networks from the segmented 3D images. The pores and throats are represented by balls and cylinders, respectively. Figure 6 illustrates the pore network structures of the three samples; in particular, the large red balls represent some of the large pores of each sample. For instance, the large red balls of the SD sample feature size between 10 and 20 μm are also shown and supported by the pore size distribution of Figure 7a and are indicated in the mean pore size of Table 3. Meanwhile, those of the A-4 and TC samples present small pore sizes between 2 and 3 μm, as shown in Figure 7b,c and Table 3.
The results of the PNM analysis are summarized in Table 3. The accuracy of the PNM analysis is depicted in Table 4, where the estimated porosities and permeabilities are compared with those obtained through laboratory experiments. The relative errors between the estimated and measured values are satisfactory and range from 3% to 10%.  1 The average number of throats connected to a pore. 2 The ratio of pore size to the size of a linked throat.

Pore Network Model
We used the pore-based algorithm in the PerGeos software to extract pore networks from the segmented 3D images. The pores and throats are represented by balls and cylinders, respectively. Figure 6 illustrates the pore network structures of the three samples; in particular, the large red balls represent some of the large pores of each sample. For instance, the large red balls of the SD sample feature size between 10 and 20 µm are also shown and supported by the pore size distribution of Figure 7a and are indicated in the mean pore size of Table 3. Meanwhile, those of the A-4 and TC samples present small pore sizes between 2 and 3 µm, as shown in Figure 7b,c and Table 3.
The results of the PNM analysis are summarized in Table 3. The accuracy of the PNM analysis is depicted in Table 4, where the estimated porosities and permeabilities are compared with those obtained through laboratory experiments. The relative errors between the estimated and measured values are satisfactory and range from 3% to 10%.

Pore Throat Size Analysis
The pore network model provides the pore and throat size distributions, as shown in Figure 7. This figure indicates that the throat distribution is wider than that of the pores for the three samples; this is consistent with the values of the number of throats extracted from PNM (see Table 3). Figure 7a shows that the dominant pore and throat sizes in the SD sample are approximately 12 µm and 5 µm, respectively, while Figure 7b,c indicates that the A-4 and TC samples feature similar pore and throat sizes, and their prevailing sizes are approximately 2 µm and 0.5 µm for the A-4 and TC samples, respectively. It is worth highlighting that the range of throat sizes extracted using PNM coincides with those obtained experimentally ( Figure 3) for the SD and A-4 samples; these span between 0.2 and 30 µm and 0.02 and 5 µm, respectively. For the TC sample, the range detected by PNM spans between 0.05 and 5 µm, while that obtained experimentally spans between 0.001 and 40 µm. The reason is that the 500 × 500 voxel, extracted from the central region of the raw image, was not enough to include big pores (see Figure 2c).

Topological Structure Analysis
The topological structure of the rock sample can be defined by the coordination number. This corresponds to the number of throats connected to each pore [11]. Figure 8 indicates that for the three samples, more than 98% of the total pores have coordination numbers lying between 1 and 12. Figure 8 also indicates that the pores with a coordination number of 3 are the ones that occur the most. Their frequency ranges from approximately 33% to 40%. The pores with a coordination number ranging from 1 to 7 constitute more than 90% of the total pores. The frequency of the pores with a coordination number of 1 in the networks is above 11% for the three samples.

Pore Throat Size Analysis
The pore network model provides the pore and throat size distributions, as shown in Figure 7. This figure indicates that the throat distribution is wider than that of the pores for the three samples; this is consistent with the values of the number of throats extracted from PNM (see Table 3). Figure 7a shows that the dominant pore and throat sizes in the SD sample are approximately 12 μm and 5 μm, respectively, while Figure 7b,c indicates that the A-4 and TC samples feature similar pore and throat sizes, and their prevailing sizes are approximately 2 μm and 0.5 μm for the A-4 and TC samples, respectively. It is worth highlighting that the range of throat sizes extracted using PNM coincides with those obtained experimentally (Figure 3) for the SD and A-4 samples; these span between 0.2 and 30 μm and 0.02 and 5 μm, respectively. For the TC sample, the range detected by PNM spans between 0.05 and 5 μm, while that obtained experimentally spans between 0.001 and 40 μm. The reason is that the 500 × 500 voxel, extracted from the central region of the raw image, was not enough to include big pores (see Figure 2c).

Topological Structure Analysis
The topological structure of the rock sample can be defined by the coordination number. This corresponds to the number of throats connected to each pore [11]. Figure 8 indicates that for the three samples, more than 98% of the total pores have coordination numbers lying between 1 and 12. Figure 8 also indicates that the pores with a coordination number of 3 are the ones that occur the most. Their frequency ranges from approximately 33% to 40%. The pores with a coordination number ranging from 1 to 7 constitute more than 90% of the total pores. The frequency of the pores with a coordination number of 1 in the networks is above 11% for the three samples.

Aspect Ratio
The aspect ratio is defined as the ratio of pore size to the size of a linked throat. The aspect ratios of the three rock samples are shown in Figure 9. This figure indicates that in the three samples, only 20% of the pores and throats have the same sizes, and most of the pores are larger than the throats. In addition, Table 3 indicates that the SD sample has a larger mean average aspect ratio compared to the A-4 and TC samples. A small aspect ratio indicates that the fluid can flow through the rock structure with relatively lower resistance.

Aspect Ratio
The aspect ratio is defined as the ratio of pore size to the size of a linked throat. The aspect ratios of the three rock samples are shown in Figure 9. This figure indicates that in the three samples, only 20% of the pores and throats have the same sizes, and most of the pores are larger than the throats. In addition, Table 3 indicates that the SD sample has a larger mean average aspect ratio compared to the A-4 and TC samples. A small aspect ratio indicates that the fluid can flow through the rock structure with relatively lower resistance.

Pore Throat Shape Factor
Accurate characterizations of pore and throat geometries are necessary for the simulation of two-phase flow in porous media. Since the shapes of real pore and throat crosssections are generally highly complex and irregular, it is essential to simplify the pores and throats into geometric bodies with idealized cross-sectional shapes. As shown in Figure 10, real pores are replaced by circles, triangles, and square pores in most pore network algorithms [35,36]. The shape factor is introduced to match the real pore space structure, which is defined as G = A/P 2 , where A and P indicate the cross-section area and the perimeter of the pore space, respectively. The shape factor identification principles are as follows: A larger shape factor value means a smoother structure. Figure 11 shows the pore and throat shape factor distribution. The figure indicates that there are a few pores and throats with G > 1/16. Hence, these can be replaced by cylinders (circles in 2D). The figure also indicates that a shape factor of approximately 12-31% of pores and 14-20% of throats is within the interval [√3/36, 1/16]. Therefore, they can be replaced by squares. Around 69-85% of pores and 76-84% of throats can be replaced by triangles.   Figure 9. Aspect ratio distribution frequency.

Pore Throat Shape Factor
Accurate characterizations of pore and throat geometries are necessary for the simulation of two-phase flow in porous media. Since the shapes of real pore and throat cross-sections are generally highly complex and irregular, it is essential to simplify the pores and throats into geometric bodies with idealized cross-sectional shapes. As shown in Figure 10, real pores are replaced by circles, triangles, and square pores in most pore network algorithms [35,36]. The shape factor is introduced to match the real pore space structure, which is defined as G = A/P 2 , where A and P indicate the cross-section area and the perimeter of the pore space, respectively. The shape factor identification principles are as follows: − Circle i f G > 1/16 A larger shape factor value means a smoother structure. Figure 11 shows the pore and throat shape factor distribution. The figure indicates that there are a few pores and throats with G > 1/16. Hence, these can be replaced by cylinders (circles in 2D). The figure also indicates that a shape factor of approximately 12-31% of pores and 14-20% of throats is within the interval √ 3/36, 1/16 . Therefore, they can be replaced by squares. Around 69-85% of pores and 76-84% of throats can be replaced by triangles.

Pore Throat Shape Factor
Accurate characterizations of pore and throat geometries are necessary for the simulation of two-phase flow in porous media. Since the shapes of real pore and throat crosssections are generally highly complex and irregular, it is essential to simplify the pores and throats into geometric bodies with idealized cross-sectional shapes. As shown in Figure 10, real pores are replaced by circles, triangles, and square pores in most pore network algorithms [35,36]. The shape factor is introduced to match the real pore space structure, which is defined as G = A/P 2 , where A and P indicate the cross-section area and the perimeter of the pore space, respectively. The shape factor identification principles are as follows: A larger shape factor value means a smoother structure. Figure 11 shows the pore and throat shape factor distribution. The figure indicates that there are a few pores and throats with G > 1/16. Hence, these can be replaced by cylinders (circles in 2D). The figure also indicates that a shape factor of approximately 12-31% of pores and 14-20% of throats is within the interval [√3/36, 1/16]. Therefore, they can be replaced by squares. Around 69-85% of pores and 76-84% of throats can be replaced by triangles.

Pore Throat Shape Factor
Accurate characterizations of pore and throat geometries are necessary for the simulation of two-phase flow in porous media. Since the shapes of real pore and throat crosssections are generally highly complex and irregular, it is essential to simplify the pores and throats into geometric bodies with idealized cross-sectional shapes. As shown in Figure 10, real pores are replaced by circles, triangles, and square pores in most pore network algorithms [35,36]. The shape factor is introduced to match the real pore space structure, which is defined as G = A/P 2 , where A and P indicate the cross-section area and the perimeter of the pore space, respectively. The shape factor identification principles are as follows: A larger shape factor value means a smoother structure. Figure 11 shows the pore and throat shape factor distribution. The figure indicates that there are a few pores and throats with G > 1/16. Hence, these can be replaced by cylinders (circles in 2D). The figure also indicates that a shape factor of approximately 12-31% of pores and 14-20% of throats is within the interval [√3/36, 1/16]. Therefore, they can be replaced by squares. Around 69-85% of pores and 76-84% of throats can be replaced by triangles.   Figure 11. (a) Pore and (b) throat shape factor distribution function.

Primary Drainage
The capillary pressure, needed to validate the flow simulations of the water-oil system, is given as MICP curves in Figure 3. Therefore, these curves, obtained from mercury-air primary drainage experiments, need to be converted to the capillary pressure in the oilwater system as follows: The capillary pressure is defined as: where σ is the interfacial tension, θ is the contact angle, defined as the angle between the liquid-vapor and solid-liquid interfaces of a solid-liquid-vapor system, and r is the throat radius. Pressure conversion is conducted using the following equation: where the subscript ow represents the oil and water system, while ha refers to the mercury and air system. The data used in the experiments and simulations are summarized in Table 5. An oil-water system was considered flowing through the pore network models. We assumed an oil-water interfacial tension of 32 Dynes/cm and a receding contact angle of 30 • . The primary drainage capillary pressures were calculated and compared with the converted MICP results of Figure 3, as shown in Figure 12. For the SD sample, the calculated capillary pressure curve matches the MICP curve. This is because the used micro-CT image resolution is sufficient to resolve all of the pores and throats of the SD sample, resulting in the correct prediction of the simulated capillary pressure. However, for A-4 and TC, the simulated capillary pressure curves are below the MICP curves at low water saturation. This is because the image resolution was not enough to capture the nanopores sized below the 0.8 micro-meter image resolution, while some of these pores and throats contribute to the experimentally determined capillary pressure values. Indeed, Figure 12, which depicts comparisons between laboratory and simulated capillary pressures obtained from high-and low-resolution images, indicates that resolving smaller pores results in a better estimation of the MICP curves.

Primary Drainage
The capillary pressure, needed to validate the flow simulations of the water-oil system, is given as MICP curves in Figure 3. Therefore, these curves, obtained from mercuryair primary drainage experiments, need to be converted to the capillary pressure in the oil-water system as follows: The capillary pressure is defined as: where σ is the interfacial tension, θ is the contact angle, defined as the angle between the liquid-vapor and solid-liquid interfaces of a solid-liquid-vapor system, and r is the throat radius. Pressure conversion is conducted using the following equation: = where the subscript ow represents the oil and water system, while ha refers to the mercury and air system. The data used in the experiments and simulations are summarized in Table 5. An oil-water system was considered flowing through the pore network models. We assumed an oil-water interfacial tension of 32 Dynes/cm and a receding contact angle of 30°. The primary drainage capillary pressures were calculated and compared with the converted MICP results of Figure 3, as shown in Figure 12. For the SD sample, the calculated capillary pressure curve matches the MICP curve. This is because the used micro-CT image resolution is sufficient to resolve all of the pores and throats of the SD sample, resulting in the correct prediction of the simulated capillary pressure. However, for A-4 and TC, the simulated capillary pressure curves are below the MICP curves at low water saturation. This is because the image resolution was not enough to capture the nanopores sized below the 0.8 micro-meter image resolution, while some of these pores and throats contribute to the experimentally determined capillary pressure values. Indeed, Figure 12, which depicts comparisons between laboratory and simulated capillary pressures obtained from high-and low-resolution images, indicates that resolving smaller pores results in a better estimation of the MICP curves. In the following, we explore the rock and two-phase flow properties using the pore network model used above to estimate the porosity, absolute permeability, capillary pressure, and relative permeabilities. Unfortunately, we do not have experimental measurements of the relative permeabilities to validate our simulated results. However, as shown below, the predicted results for relative permeability trends and values are consistent with the expected multiphase flow behavior in porous media.
For instance, the water-and oil-predicted relative permeabilities for the primary drainage with two different image resolutions (one that resolves all of the pores and throats, while the other does not) are presented in Figure 13a-c for SD, A-4, TC, respectively. The low image resolution indicates that the oil starts to flow at a water saturation of around 0.85, while the higher-resolution one shows that oil starts to flow at a water saturation of around 0.65. Since the lower resolution captures only relatively larger-diameter throats, the injected oil encounters less resistance to flow in pore throats and thus starts to invade the water at higher saturations. However, high image resolution allows capturing smaller pore throats, the water-being the wetting phase-preferentially occupies these throats. This induces a higher resistance to the oil flow in the rock sample, resulting in a delay of oil movement in the rock until a 0.65 water saturation. Overall, Figure  13 shows the oil relative permeability increases with a higher image resolution. This is because a higher image resolution can detect more throat branches connected to the resolved and smaller pores, resulting in higher relative permeability values.
Furthermore, the results in Figure 13 from higher resolution images indicate that the water relative permeability drops significantly with the oil injection, while the oil relative permeability remains close to its initial zero value. This can be explained by the fact that the water phase occupies most of the pores and throats in the network model. With the increase in oil saturation (decrease in water saturation), the oil relative permeability increases gradually. The water relative permeability values reached almost 0 when the oil saturation exceeded 0.5. The very low but non-zero water relative permeability values indicate that the injected oil did not completely block the water pathways. The water continues to flow through the connected wetting water layers along the surfaces of the pore space. At the cross-over saturation, both the oil and water relative permeability values are around 0.1 or lower. Such low relative permeabilities are explained by the presence of arc menisci, which act as flow barriers to the flow of the two phases [37]. In the following, we explore the rock and two-phase flow properties using the pore network model used above to estimate the porosity, absolute permeability, capillary pressure, and relative permeabilities. Unfortunately, we do not have experimental measurements of the relative permeabilities to validate our simulated results. However, as shown below, the predicted results for relative permeability trends and values are consistent with the expected multiphase flow behavior in porous media.
For instance, the water-and oil-predicted relative permeabilities for the primary drainage with two different image resolutions (one that resolves all of the pores and throats, while the other does not) are presented in Figure 13a-c for SD, A-4, TC, respectively. The low image resolution indicates that the oil starts to flow at a water saturation of around 0.85, while the higher-resolution one shows that oil starts to flow at a water saturation of around 0.65. Since the lower resolution captures only relatively larger-diameter throats, the injected oil encounters less resistance to flow in pore throats and thus starts to invade the water at higher saturations. However, high image resolution allows capturing smaller pore throats, the water-being the wetting phase-preferentially occupies these throats. This induces a higher resistance to the oil flow in the rock sample, resulting in a delay of oil movement in the rock until a 0.65 water saturation. Overall, Figure 13 shows the oil relative permeability increases with a higher image resolution. This is because a higher image resolution can detect more throat branches connected to the resolved and smaller pores, resulting in higher relative permeability values.
Furthermore, the results in Figure 13 from higher resolution images indicate that the water relative permeability drops significantly with the oil injection, while the oil relative permeability remains close to its initial zero value. This can be explained by the fact that the water phase occupies most of the pores and throats in the network model. With the increase in oil saturation (decrease in water saturation), the oil relative permeability increases gradually. The water relative permeability values reached almost 0 when the oil saturation exceeded 0.5. The very low but non-zero water relative permeability values indicate that the injected oil did not completely block the water pathways. The water continues to flow through the connected wetting water layers along the surfaces of the pore space. At the cross-over saturation, both the oil and water relative permeability values are around 0.1 or lower. Such low relative permeabilities are explained by the presence of arc menisci, which act as flow barriers to the flow of the two phases [37].

Waterflooding
After the primary drainage, we conducted simulations of water injection into the pore networks (water flooding). Since it is widely accepted that carbonate rocks are oilwater systems, we assumed that the samples are slightly oil-wet. This assumption is reasonable because oil migrates into the pore space and displaces water in the rocks, which are fully saturated with water at their formation. In other words, long-time exposure to oil alters the rock's wettability; the reservoir becomes an oil-water system. For the water flooding simulation, we set the advancing contact angles, which are the maximum contact angles produced during a wetting process, from 140° to 150°. We also assumed that the wettability distribution is independent of the pore size. The predicted water flooding oilwater relative permeabilities of all three samples are shown in Figure 14.
In the oil-water system, one can notice a rapid increase in water relative permeability above a certain water saturation. In the case of the SD sample, the relative permeability increased rapidly when the water saturation reached 0.4. This indicates a rapid water connection and higher water flow conductance. At the end of the water flooding process, the water relative permeabilities reached almost 1 for the three samples. This means that the water phase (non-wetting phase) occupies the center of pore space, while the wetting oil is stranded in small pore spaces. Another result is that the oil relative permeability dropped to values lower than 0.01 after a water saturation of 0.5 for the three samples. The exact values are 51.32%, 67.83%, and 53.24% for the SD, A-4, and TC samples, respectively. The oil relative permeability remaining at very low values until the end of waterflooding is the signature of layer drainage [30,37]. This means that the connected layers of the oil phase with very low conductivity prevent the trapping of the oil in the pore space.

Waterflooding
After the primary drainage, we conducted simulations of water injection into the pore networks (water flooding). Since it is widely accepted that carbonate rocks are oil-water systems, we assumed that the samples are slightly oil-wet. This assumption is reasonable because oil migrates into the pore space and displaces water in the rocks, which are fully saturated with water at their formation. In other words, long-time exposure to oil alters the rock's wettability; the reservoir becomes an oil-water system. For the water flooding simulation, we set the advancing contact angles, which are the maximum contact angles produced during a wetting process, from 140 • to 150 • . We also assumed that the wettability distribution is independent of the pore size. The predicted water flooding oil-water relative permeabilities of all three samples are shown in Figure 14.
In the oil-water system, one can notice a rapid increase in water relative permeability above a certain water saturation. In the case of the SD sample, the relative permeability increased rapidly when the water saturation reached 0.4. This indicates a rapid water connection and higher water flow conductance. At the end of the water flooding process, the water relative permeabilities reached almost 1 for the three samples. This means that the water phase (non-wetting phase) occupies the center of pore space, while the wetting oil is stranded in small pore spaces. Another result is that the oil relative permeability dropped to values lower than 0.01 after a water saturation of 0.5 for the three samples. The exact values are 51.32%, 67.83%, and 53.24% for the SD, A-4, and TC samples, respectively. The oil relative permeability remaining at very low values until the end of waterflooding is the signature of layer drainage [30,37]. This means that the connected layers of the oil phase with very low conductivity prevent the trapping of the oil in the pore space.

The Effect of Wettability on Relative Permeability
In this section, we present the results of our investigations of the effect of rock wettability on the relative permeability in two different wettability systems, namely, uniform and mixed wet cases. In the uniform wet system, the same wettability is assigned to all of the pores and throats of the rock sample pore network [32]. Meanwhile, in the mixed wet system, some pores and throats are considered water-wet, while others are oil-wet. We considered water flooding of an initially oil-saturated rock sample, and wettability is selected by using the appropriate maximum and minimum advancing contact angle values with a difference of 10° within the PerGeos pore network model, knowing that other studies have considered a difference of 20° [32]. We used a difference of 20° and the obtained variation of relative permeability between different wettability conditions was not easily captured. While with the selected 10° difference, this variation is observed as will be shown in the relative permeability curves below.
The literature indicates that in the water-wet system the advancing contact angle is less than 90 degrees while in the oil-water system the advancing contact angle is more than 140 degrees [38]. We used the Amott-Harvey index to determine the system's wettability. The Amott-Harvey index is the difference between the Amott water index and the Amott oil index [38]. The Amott water (oil) index is defined as the ratio of the amount of oil (water) displaced by spontaneous imbibition of water (oil) to the total amount of oil (water) displaced during imbibition (or secondary drainage) and is calculated by: where Sspp is the phase p saturation displaced by spontaneous imbibition, Swi is the initial water saturation, and Sor is the residual oil saturation.

Uniform Wet System
To investigate the effect of wettability on the Amott-Harvey index, we varied the advancing contact angle from 40° to 165°. Within the PerGeo package and specifically the

The Effect of Wettability on Relative Permeability
In this section, we present the results of our investigations of the effect of rock wettability on the relative permeability in two different wettability systems, namely, uniform and mixed wet cases. In the uniform wet system, the same wettability is assigned to all of the pores and throats of the rock sample pore network [32]. Meanwhile, in the mixed wet system, some pores and throats are considered water-wet, while others are oil-wet. We considered water flooding of an initially oil-saturated rock sample, and wettability is selected by using the appropriate maximum and minimum advancing contact angle values with a difference of 10 • within the PerGeos pore network model, knowing that other studies have considered a difference of 20 • [32]. We used a difference of 20 • and the obtained variation of relative permeability between different wettability conditions was not easily captured. While with the selected 10 • difference, this variation is observed as will be shown in the relative permeability curves below.
The literature indicates that in the water-wet system the advancing contact angle is less than 90 degrees while in the oil-water system the advancing contact angle is more than 140 degrees [38]. We used the Amott-Harvey index to determine the system's wettability. The Amott-Harvey index is the difference between the Amott water index and the Amott oil index [38]. The Amott water (oil) index is defined as the ratio of the amount of oil (water) displaced by spontaneous imbibition of water (oil) to the total amount of oil (water) displaced during imbibition (or secondary drainage) and is calculated by: where S spp is the phase p saturation displaced by spontaneous imbibition, S wi is the initial water saturation, and S or is the residual oil saturation.

Uniform Wet System
To investigate the effect of wettability on the Amott-Harvey index, we varied the advancing contact angle from 40 • to 165 • . Within the PerGeo package and specifically the pore network model, to select an advancing contact angle, a minimum value and maximum value have to be specified rather than a constant value of the angle. In their work on wettability effects on water flooding at the pore scale, Zhao et al. [32] used a difference of 20 • between these minimum and maximum values. We carried out some simulations with both 20 • and 10 • differences and found no pronounced variation in the obtained results.
Therefore, we selected a difference of 10 • to stay closer to the specific advancing contact angle value in all of our simulations. The computed Amott-Harvey index is presented in Figure 15 below. This figure shows that the index decreases from 1 (water-wet system) to −1 (oil-water system) within the used advancing contact angle range from 40 • to 165 • . Furthermore, the Amott-Harvey index went through an intermediate plateau of zero between the advancing contact angle of the water-wet system (90 • ) and the oil-water system (140 • ), indicating natural wettability behavior in this range [38]. pore network model, to select an advancing contact angle, a minimum value and maximum value have to be specified rather than a constant value of the angle. In their work on wettability effects on water flooding at the pore scale, Zhao et al. [32] used a difference of 20° between these minimum and maximum values. We carried out some simulations with both 20° and 10° differences and found no pronounced variation in the obtained results. Therefore, we selected a difference of 10° to stay closer to the specific advancing contact angle value in all of our simulations. The computed Amott-Harvey index is presented in Figure 15 below. This figure shows that the index decreases from 1 (water-wet system) to −1 (oil-water system) within the used advancing contact angle range from 40° to 165°. Furthermore, the Amott-Harvey index went through an intermediate plateau of zero between the advancing contact angle of the water-wet system (90°) and the oil-water system (140°), indicating natural wettability behavior in this range [38]. Next, we investigated the effect of the advancing contact angle (wettability) on the relative permeability curves in the uniform water-wet and oil-water systems, respectively, during water flooding (initially fully oil-saturated rock sample). The predicted values of the oil-water relative permeabilities of the three samples for different advancing contact angles in the water-wet and oil-water systems are shown in the right and left columns of Figure 16, respectively. In each of these uniform wet systems, the wetting phase is known to occupy the smaller pores and throats and to form thin films attached to the walls of the larger pores and throats. On the contrary, the non-wetting phase fills the central volume of the pores. This last behavior influences the interpretation of the relative permeability curves.
The left column in Figure 16 shows the influence of the wettability on the relative permeabilities in a water-wet system (advancing contact angle less than 90). Figure 16a for the SD sample shows that increasing the advancing contact angle in a water-wet system has almost no influence on both the water and oil relative permeability curves. In addition, the oil relative permeability sharply decreased as soon as the injected water saturation started increasing in the system, reaching a value of almost zero at a water saturation of approximately 0.3. Meanwhile, the water relative permeability remained almost zero until a water saturation of 0.2, then increased to a value of only 0.025 at a water saturation of 0.3, and both relative permeabilities remain unchanged with increasing water injection (water saturation). Beyond a water saturation of 0.3, the oil was trapped in the central volume of pores and the further injected water slipped over it by adhering to the water-wet pore and throat surfaces. This was also a consequence of the fact that most of the pores in the SD sample were poorly connected. For instance, 27% of the total pores in this sample and 60% of the pores with a radius smaller than 10 μm were only connected to one pore, as shown in Figures 8 and 17, respectively. Figure 16c,e correspond to the A-4 and TC samples, respectively, both showing that the oil relative permeability decreased with increasing advancing contact angle. They also show that the oil relative permeability decreased with water saturation, but not as sharply Next, we investigated the effect of the advancing contact angle (wettability) on the relative permeability curves in the uniform water-wet and oil-water systems, respectively, during water flooding (initially fully oil-saturated rock sample). The predicted values of the oil-water relative permeabilities of the three samples for different advancing contact angles in the water-wet and oil-water systems are shown in the right and left columns of Figure 16, respectively. In each of these uniform wet systems, the wetting phase is known to occupy the smaller pores and throats and to form thin films attached to the walls of the larger pores and throats. On the contrary, the non-wetting phase fills the central volume of the pores. This last behavior influences the interpretation of the relative permeability curves.
The left column in Figure 16 shows the influence of the wettability on the relative permeabilities in a water-wet system (advancing contact angle less than 90). Figure 16a for the SD sample shows that increasing the advancing contact angle in a water-wet system has almost no influence on both the water and oil relative permeability curves. In addition, the oil relative permeability sharply decreased as soon as the injected water saturation started increasing in the system, reaching a value of almost zero at a water saturation of approximately 0.3. Meanwhile, the water relative permeability remained almost zero until a water saturation of 0.2, then increased to a value of only 0.025 at a water saturation of 0.3, and both relative permeabilities remain unchanged with increasing water injection (water saturation). Beyond a water saturation of 0.3, the oil was trapped in the central volume of pores and the further injected water slipped over it by adhering to the water-wet pore and throat surfaces. This was also a consequence of the fact that most of the pores in the SD sample were poorly connected. For instance, 27% of the total pores in this sample and 60% of the pores with a radius smaller than 10 µm were only connected to one pore, as shown in Figures 8 and 17, respectively. Figure 16c,e correspond to the A-4 and TC samples, respectively, both showing that the oil relative permeability decreased with increasing advancing contact angle. They also show that the oil relative permeability decreased with water saturation, but not as sharply as that of the SD sample and reached a value of approximately zero at water saturation of approximately 0.5 compared to the 0.3 value of the SD sample. change with the increase in advancing contact angle. The explanation is that in the oilwater system, the connection of oil layers has already been established. The path of the oil flow through the medium cannot be blocked by the injected water. It was noticed that the oil relative permeability reached almost zero at a water saturation 0.5 and remained at this value with the increasing water injection. This means that the connected oil layers had very low conductance at high water saturation, hence contributing to residual oil saturation.

Mixed Wet System
The effect of wettability on the relative permeabilities in a mixed wet system was studied by varying the fraction of oil-wet pores from 0.2 to 0.8, while the advancing contact angle was fixed between a minimum and maximum value of 30° and 40° for water and of 160° and 170° for oil, during waterflooding of the initially fully oil-saturated sample. The relative permeability curves are presented in Figure 18 for the three samples. Note that the fraction with a value 0 (1) corresponds to a completely water-wet system (oil-wet). In addition, Figure 16c shows that increasing the advancing contact angle leads to a decrease in the water relative permeability. This can be explained by the presence of poorly connected water paths; 50% of the pores smaller than 10 µm had one, two, and three throat connections, which also resulted in a sharp decrease in the oil relative permeability.
Meanwhile, for the TC sample and opposite to the A-4 sample, Figure 16e shows that the water relative permeability increased with the increase in the advancing contact angle. This is probably due to the fact that the pore and throat sizes in the former were almost double those of the latter, and the absolute permeability was much higher for TC than for A-4 (see Table 2). Thus, the water preferentially invaded larger pores/throats with better conductance with decreasing water wettability.
In the oil-water system (right column of Figure 16), with advancing contact angles ranging from 140 to 175, the relative permeabilities of both the oil and water did not change with the increase in advancing contact angle. The explanation is that in the oil-water system, the connection of oil layers has already been established. The path of the oil flow through the medium cannot be blocked by the injected water. It was noticed that the oil relative permeability reached almost zero at a water saturation 0.5 and remained at this value with the increasing water injection. This means that the connected oil layers had very low conductance at high water saturation, hence contributing to residual oil saturation.

Mixed Wet System
The effect of wettability on the relative permeabilities in a mixed wet system was studied by varying the fraction of oil-wet pores from 0.2 to 0.8, while the advancing contact angle was fixed between a minimum and maximum value of 30 • and 40 • for water and of 160 • and 170 • for oil, during waterflooding of the initially fully oil-saturated sample. The relative permeability curves are presented in Figure 18 for the three samples. Note that the fraction with a value 0 (1) corresponds to a completely water-wet system (oil-wet). The fraction of oil-wet pores and throats within PerGeos was assigned by randomly selecting the corresponding number of pores and throats and setting their respective advancing contact angle to their water-or oil-wet value. Figure 18 shows that the oil relative permeability consistently decreased with an increasing fraction of oil-wet pores and throats (left side of Figure 18a-c, corresponding to lower water saturation). Meanwhile, we noticed the opposite trend in the water relative permeability (right side od Figure 18a-c, corresponding to higher water saturation). On the lower water saturation side of the three figures, we noticed a sharper decrease of the oil relative permeability for the SD sample than that for A-4 and TC. This resulted in the oil relative permeability becoming almost zero at a water saturation of approximately 0.2 for the SD sample, while reaching this value at a water saturation of approximately 0.4 for the A-4 and TC samples. The last observations are in agreement with the earlier discussed water-wet and oil-water systems shown in Figure 16 above.

Mixed Wet System
The effect of wettability on the relative permeabilities in a mixed wet system was studied by varying the fraction of oil-wet pores from 0.2 to 0.8, while the advancing contact angle was fixed between a minimum and maximum value of 30° and 40° for water and of 160° and 170° for oil, during waterflooding of the initially fully oil-saturated sample. The relative permeability curves are presented in Figure 18 for the three samples. Note that the fraction with a value 0 (1) corresponds to a completely water-wet system (oil-wet). The fraction of oil-wet pores and throats within PerGeos was assigned by randomly selecting the corresponding number of pores and throats and setting their respective advancing contact angle to their water-or oil-wet value. Figure 18 shows that the oil relative permeability consistently decreased with an increasing fraction of oil-wet pores and throats (left side of Figure 18a-c, corresponding to lower water saturation). Meanwhile, we noticed the opposite trend in the water relative permeability (right side od Figure 18a-c, corresponding to higher water saturation). On the lower water saturation side of the three figures, we noticed a sharper decrease of the oil relative permeability for the SD sample than that for A-4 and TC. This resulted in the oil relative permeability becoming almost zero at a water saturation of approximately 0.2 for the SD sample, while reaching this value at a water saturation of approximately 0.4 for the A-4 and TC samples. The last observations are in agreement with the earlier discussed water-wet and oil-water systems shown in Figure 16 above.

Conclusions
In this study, we applied the pore network model to investigate carbonate rock properties in terms of structure, morphology, and multiphase flow properties. In the beginning, we presented the experimental micro-CT images and laboratory measurement results from three different carbonate rock samples that featured quite different pore and throat distributions, capillary pressure, and normalized permeability. With the appropriate image resolution, we were able to predict values for the network-connected porosity, total porosity, and absolute permeability that agree well with those measured experimentally. In addition, when the image resolution was good enough to resolve the relevant pore and throat sizes, the simulated capillary pressure curves matched well with the measured MICP curves. It is worth highlighting that the simulated capillary pressure curves were obtained without matching any parameters such as contact angle or pore size because of the use of the appropriate resolution of the rock images. The relative permeability curves were substantially affected by the micro-CT image resolution. Higher image resolution allowed to capture smaller pores and throats, which control the multiphase flow phenomena in porous media. We also studied the effect of wettability on the relative permeabilities both in water-wet and oil-water systems through the variation in the advancing contact angle, and we observed a pronounced effect of wettability on relative permeability curves for all three investigated carbonate rock samples. In addition, in the mixed wet system, we investigated the effect of the fraction of the oil-wet pores and throats, and we also noted a strong influence of this parameter on the relative permeability curves.   Figure 18. Effect of the oil-wet fraction on the oil-water relative permeability in the mixed wet system of (a) SD, (b) A-4, and (c) TC.

Conclusions
In this study, we applied the pore network model to investigate carbonate rock properties in terms of structure, morphology, and multiphase flow properties. In the beginning, we presented the experimental micro-CT images and laboratory measurement results from three different carbonate rock samples that featured quite different pore and throat distributions, capillary pressure, and normalized permeability. With the appropriate image resolution, we were able to predict values for the network-connected porosity, total porosity, and absolute permeability that agree well with those measured experimentally. In addition, when the image resolution was good enough to resolve the relevant pore and throat sizes, the simulated capillary pressure curves matched well with the measured MICP curves. It is worth highlighting that the simulated capillary pressure curves were obtained without matching any parameters such as contact angle or pore size because of the use of the appropriate resolution of the rock images. The relative permeability curves were substantially affected by the micro-CT image resolution. Higher image resolution allowed to capture smaller pores and throats, which control the multiphase flow phenomena in porous media. We also studied the effect of wettability on the relative permeabilities both in water-wet and oil-water systems through the variation in the advancing contact angle, and we observed a pronounced effect of wettability on relative permeability curves for all three investigated carbonate rock samples. In addition, in the mixed wet system, we investigated the effect of the fraction of the oil-wet pores and throats, and we also noted a strong influence of this parameter on the relative permeability curves.