The wettability of mineral surfaces has significant impacts on transport mechanisms of two-phase flow, distribution characteristics of fluids, and the formation mechanisms of residual oil during water flooding. However, few studies have investigated such effects of mineral type and its surface wettability on rock properties in the literature. To unravel the dependence of hydrodynamics on wettability and minerals distribution, we designed a new experimental procedure that combined the multiphase flow experiments with a CT scan and QEMSCAN to obtain 3D digital models with multiple minerals and fluids. With the aid of QEMSCAN, six mineral components and two fluids in sandstones were segmented from the CT data based on the histogram threshold and watershed methods. Then, a mineral surface analysis algorithm was proposed to extract the mineral surface and classify its mineral categories. The in situ contact angle and pore occupancy were calculated to reveal the wettability variation of mineral surface and distribution characteristics of fluids. According to the shape features of the oil phase, the self-organizing map (SOM) method, one of the machine learning methods, was used to classify the residual oil into five types, namely, network, cluster, film, isolated, and droplet oil. The results indicate that each mineral’s contribution to the mineral surface is not proportional to its relative content. Feldspar, quartz, and clay are the main minerals in the studied sandstones and play a controlling role in the wettability variation. Different wettability samples show various characteristics of pore occupancy. The water flooding front of the weakly water-wet to intermediate-wet sample is uniform, and oil is effectively displaced in all pores with a long oil production period. The water-wet sample demonstrates severe fingering, with a high pore occupancy change rate in large pores and a short oil production period. The residual oil patterns gradually evolve from networks to clusters, isolated, and films due to the effects of snap-off and wettability inversion. This paper reveals the effects of wettability of mineral surface on the distribution characteristics and formation mechanisms of residual oil, which offers us an in-deep understanding of the impacts of wettability and minerals on multiphase flow and helps us make good schemes to improve oil recovery.

Enhancing oil recovery (EOR) is crucial for oilfields to obtain stable oil production. Water flooding is a dominating enhanced oil recovery technique owing to its low cost, easy operation, and less pollution [14]. Therefore, it is significant for EOR to study the transport mechanisms of multiphase flow and distribution characteristics of fluids. Pore-scale fluid flow experiment is an effective method to investigate the transport mechanisms and distribution characteristics of residual oil, including microfluidics [58], natural sandstone models [913], nuclear magnetic resonance [1416], confocal laser [1719], and CT scan [2022]. Owing to the rapid development of CT technology, it has been widely used in EOR [2327], carbon dioxide storage [28], environmental governance [29], water-rock interaction [30, 31], reservoir modelling [3236], hydraulic conductivity [37, 38], and reservoir evaluation [3947]. Due to the advantages of high resolution, nondestructive, and in situ, CT scan has become an effective method for investigating residual oil. For example, Wang et al. carried out displacement experiments of water flooding, polymer flooding, preformed particle gel (PPG) flooding, and obtained CT images of different displacement stages [48]. Then the content, distribution characteristics, and morphological characteristics of residual oil under different displacement systems were analyzed. An et al. designed in situ displacement experiments based on CT scans and reveal residual oil distribution after different displacement methods and stages [49]. The shape of the residual oil was further classified according to shape factor and Euler number. Guo et al. visually and quantitatively revealed the effects of pore structure and permeability on oil displacement efficiency and residual oil morphology in water-wet sandstone samples with high-resolution micro-computer tomography (μCT) technology [50]. The fractal dimension method was used to quantitatively investigate the complexity of the pore structure of each sample and the distribution mechanism of the residual oil. Wu et al. carried out the two-phase (brine and oil) flow experiment on a water-wet sandstone based on CT scans. A comprehensive analysis of the non-wetting (oil) phase distribution in 3D pore space and individual pores is presented [20]. The axis-aligned bounding box (AABB) algorithm was proposed to determine which oil blobs exist in a unique pore. Lin et al. [51] proposed core flooding experiments integrated with X-ray microtomography (micro-CT) to visualize asphaltene precipitation in the pore space of rocks and assess the reduction in permeability. The results displayed that the asphaltene tended to clog the larger pores, causing a 29-fold reduction in permeability. Qin et al. investigated near-miscible supercritical CO2 (scCO2) injection and its underlying pore-scale physics in an oil-wet carbonate rock at elevated temperature and pressure conditions using micro-CT [23]. The results indicated that the scCO2 injection triggered a wettability reversal process resulting in reduced oil wetness. Consequently, pore sizes neither dictated any preferential invasion order nor restricted the displacement efficiency of the injection process. Lin et al. applied 3D micro-CT to image a capillary drainage process in a cm-scale heterogeneous laminated sandstone containing three distinct regions with different pore sizes to study the capillary pressure [52]. They illustrated that fluid displacement is controlled by pore size, connectivity, and proximity to the boundaries between regions with different pore-space characteristics. Ju et al. used X-ray computed tomography, integrated with triaxial loading techniques, to capture in situ the immiscible water-oil displacement and oil trapping inside 3D pore spaces during the deformation induced by various geostresses [53]. The results confirmed that stress-induced pore structure deformation has an evident impact on the 3D water-oil displacement behavior and efficiency.

Wettability is a tendency of a fluid to spread on or adhere to a solid and mineral surface in the presence of other immiscible fluids, which is perhaps the most critical factor that affects the oil recovery [54, 55]. Wettability has a crucial impact on the behavior of multiphase flow during EOR [56]. For example, Iglauer et al. [57] analyzed the size distribution, morphology, and saturation of residual oil under different wettability conditions by X-ray microtomography. The results demonstrated that the residual oil had an apparent flat geometric structure in oil-wet rock, mainly on the particle surface and tiny pores. The residual oil was spherical in water-wet rock, primarily located in the middle of large pores. Tan et al. [58] investigated the microscopic displacement mechanism of water flooding and polymer flooding of hydrophilic and weakly lipophilic rocks by CT scan technology. The results illustrated that the flow mechanism of water flooding and polymer flooding of hydrophilic rock was mainly “crawling” along the rock’s surface, while that of weak lipophilic rock was mostly “advancing” along the middle of the pore. The most widely used wettability measurement techniques include contact angle measurement, the Amott-Harvey test, and the United States Bureau of Mines (USBM) method [59]. But these techniques cannot characterize the wettability variation of pore space [60]. The development of high-resolution imaging of three-dimensional porous media with multifluid allows the in situ measurement of contact angle [61, 62]. Based on high-resolution images, many scholars have developed a variety of methods to assess wettability, such as integral geometry [63, 64], topological principles [65], and energy balance [66, 67]. These methods can effectively evaluate the wetting properties of complex porous media. In this study, an in situ contact angle measurement method developed by AlRatrout et al. [60] was used to characterize the wettability based on X-ray imaging, which provides a rapid and accurate characterization of pore-scale wettability. In addition, the wettability of each mineral is variable, which results in wettability heterogeneities [6871]. Different mineral components show different oil-water affinities. For example, Arshadi et al. [72] indicated that quartz and feldspar grains were mostly weakly water-wet, whereas calcite and dolomite were intermediate and weakly oil-wet. Ferrari et al. [55] displayed that quartz minerals could increase the hydrophilic degree of the rock surface. But it is difficult to accurately segment different mineral components by only relying on the grayscale value of CT data since grayscale values of some minerals (e.g., quartz and feldspar) are very close.

In summary, few published studies comprehensively investigated the effects of the wettability of different mineral surfaces on multiphase flow mechanisms and oil distribution characteristics. Combining CT scan with other mineral scanning methods (such as QEMSCAN technology) is necessary to expand CT data’s mineral component information, establish the distribution law of wettability at the pore scale, and further analyze its influence on the formation and distribution of residual oil. Therefore, combined with the QEMSCAN results, the CT data was divided into multiple mineral components in this paper, and the contact angles of different mineral components were calculated to analyze wettability, and finally, the control effects of wettability on the formation mechanism and distribution pattern of the residual oil were revealed.

The rest of the paper is organized as follows (Figure 1). We introduced a new experimental design in detail in Section 2. This section also demonstrated QEMSCAN and CT data processing, and the calculation methods of rock and fluid properties, including a mineral surface segmentation algorithm, the brine in situ contact angles measurement method, pore occupancy calculation program, and the self-organizing map (SOM) technology. In Section 3, the mineral component characteristics, the variation of wettability, evolution of residual oil saturation, pore occupancy change, and residual oil pattern were investigated. Finally, some important conclusions are listed in Section 4. The paper presents a study on the wettability of mineral surfaces and its impacts on transport mechanisms of two-phase flow, distribution characteristics of fluids, and the formation mechanisms of the residual oil during imbibition, which provides a fundamental understanding of transport mechanisms of multiphase flow and makes schemes to improve oil recovery.

2.1. Samples and Reagents

Two parallel sandstone samples (8.0mm±0.01mm in diameter and 32.5mm±0.01mm in length) were used for the multiphase flow experiments under 23 MPa and 85 °C. The helium porosity of two samples were 28.2%±0.2% and 26.3%±0.2%, respectively. The nitrogen permeability of two samples measured by Darcy’s law were 387±2×103μm2 and 316±2×103μm2, respectively. The oil phase used for the experiment is the subsurface crude oil from the Jidong field after dehydrating. It has a density of 822.5 kg/m3 and a viscosity of 4.51 mPa.s, provided by the Jidong field. The brine phase is the subsurface water filtered by medium-speed filter paper. The brine density was measured to be 982.8±0.1kg/m3 at 85 °C by weighing 60 mL of the brine. The viscosity of brine was 0.595±0.02mPa.s at 85 °C. To enhance the contrast between brine and oil, 20 wt% of potassium iodide was added to the brine. The capillary number is about 1.24×106 when the brine is injected into the core.

2.2. Experimental Setup and Procedure

Figure 2 shows the experimental setup used in this study [73]. We applied carbon-fiber core clamping in this experiment. To prevent sample shaking, we replaced all tubing close to the rotary table with soft and flexible poly-ether-ether-ketone (PEEK) tubing. This kind of tub can resist high pressure and exert very little external force on the core during rotation, reducing the likelihood of sample shaking during scanning and maintaining high image accuracy. A thin polyimide heating film was utilized for core heating. This method reduced the load on the core and allowed for precise heating, resulting in a temperature control error within ±1 °C. The fluid configurations were imaged using a Versa nanoVoxel-3502E X at a voxel size of 3 μm. The X-ray energy was 120 keV, and the exposure time was 1 s. These projections were reconstructed into a 3D image using proprietary software on the Versa system. All reconstructed images were registered to the dry scan image with the same orientation. The samples’ mineralogy was determined using QEMSCAN 650 F manufactured by FEI, now Thermo Fisher Scientific, USA. QEMSCAN creates phase assemblage maps of a specimen surface scanned by a high-energy accelerated electron beam along with a predefined raster scan pattern. Low-count energy-dispersive X-ray spectra (EDX) are generated and provide information on the elemental composition at each measurement point. In combination with back-scattered electron (BSE) brightness and X-ray count rate in formation, the elemental composition is converted into mineral phases. In this study, the X-ray source was optimized at 15 kV and 1.6 nA, and a species identification protocol was used to convert the raw data with a resolution of 5 μm per pixel. The experimental procedure is as follows:

  • (i)

    The plug sample was placed in the carbon fiber holder for a CT scan of the central part of the dry core

  • (ii)

    A vacuum pump was connected to the outlet of the core holder to draw a vacuum, and then the core was injected with brine (20 wt% KI has been added) at a constant rate of 0.03 ml/min to ensure that the sample was fully saturated with brine. At the same time, the confining pressure was gradually increased to 23 MPa, the back pressure was increased to 21 MPa, and the temperature was heated to 85 °C

  • (iii)

    Crude oil was injected at a constant rate of 0.03 ml/min until no more brine was produced at the outlet for a while. Then, the crude oil was considered adequately saturated, and a CT scan was performed at the same location

  • (iv)

    Brine (with 20 wt% KI added) was injected at a constant rate of 0.03 ml/min, and the injection was stopped when displacement to 1 PV. Then, a CT scan was performed at the same location after stabilization

  • (v)

    After scanning, the brine was injected at the same rate, and the injection was stopped when the displacement reached 20 PV. Then, a CT scan was performed at the same location after stabilization

  • (vi)

    Repeat Step V, and the injection was stopped when displacement to 50 PV. Then, a CT scan was performed at the same location after stabilization

  • (vii)

    The experiment was stopped. After that, the plug sample was taken out and cut in the middle after oil washing and drying. Then, the QEMSCAN scanning was performed to obtain the mineral component information after filming and carbon plating

2.3. QEMSCAN and CT Data Processing

Based on the QEMSCAN (Quantitative Evaluation of Minerals by scanning electron microscopy) results, the ratio of each mineral pixel to the sum of all mineral pixels is calculated as the relative content of minerals. The QEMSCAN results were further used to better identify and quantify the major minerals during the segmentation of CT images. In this study, CT data segmentation progress is as follows:

  • (i)

    All three-dimensional CT data were reconstructed on PerGeos software. Based on dry scanning, CT data were aligned and checked in order to ensure that all data were in situ for easy observation and comparative analysis. 1100×1100×1100 voxels data were cropped out for further processing

  • (ii)

    Each CT data was filtered by the nonlocal median method. The CT data were segmented into solid, brine, and oil phases based on the marker-based watershed method, providing accurate segmentation results at the brine-water contact surface (demonstrated in Figure 3).

  • (iii)

    The solid label was multiplied with dry data to eliminate the pore phase, and the multiplied data was segmented into mineralogy maps by the histogram threshold method. The gray value of the 3D image obtained by CT scanning corresponds to the density of each component of the tested sample. According to this correspondence, CT data can be segmented into multimineral components. The processes were summarized as follows: (i) The gray distribution of multiplied data was statistics (Figures 4), (ii) the segmentation threshold was set according to the gray distribution, and the segmentation results were evaluated visually, and (iii) the threshold was adjusted according to the visual evaluation results and the main mineral components; this step was repeated until the segmentation results were confirmed. Finally, the solid was segmented into quartz, feldspar, dolomite, siderite, clay, and pyrite, and the segmented data was called CTsd

2.4. Calculation of Rock and Fluid Properties

2.4.1. Mineral Content and Mineral Surface Area

We developed a mineral surface segmentation algorithm to divide the surface into different mineral categories. The process is as follows: (I) The surface of all minerals was extracted from the CTsd, the x, y, and z coordinate matrix of the surface pixel points were stored, and the coordinates located at the edges of the data volume were removed. (II) According to the coordinate matrix, the mineral labels at the top, bottom, left, right, front, and back positions of the surface point were retrieved and recorded based on the CTsd. After removing the repeated labels, the mineral category matrix of the surface point was obtained, and the surface point was classified as the mineral category shown in the category matrix. If there were multicategories, the point was repeatedly classified as multimineral categories. (III) Step (II) was repeated to complete the classification of all surface points. Then, the segmented data of the mineral surface was called CTms.

2.4.2. In Situ Contact Angle Measurement

The brine in situ contact angles and the position coordinates were measured by an open-source algorithm [54, 74]. The open-source algorithm consists of the following main steps (Figure 5).
  • (i)

    Based on the segmented data, which contains the solid, brine, and oil phase, a multi-zone mesh that consists of multiple face zones was generated by the open-source algorithm. The mesh represented the interfaces of contacted phases. The set of vertices shared by all face zones represented the three-phase contact line

  • (ii)

    A volume preserving Gaussian smoothing on all vertices belonging to the generated mesh was applied, which removes the voxelized artifacts from the segmented 3D image

  • (iii)

    A local uniform curvature smoothing was applied on vertices of the fluid/fluid interface, which is consistent with capillary equilibrium

  • (iv)

    The contact angle on each vertex that belongs to the contact line set was computed using the smoothed mesh as follows:

(1)θ=πcos1n1n2,
where n1 represents a vector normal to the oil/brine interface and n2 represents a vector normal to the brine/rock interface.

A detailed description of the method can be found in the paper by AlRatrout et al. Table 1 shows the relationship between contact angle and reservoir wettability [75]. Then, we classified the mineral categories of each contact angle based on the CTms. The process is as follows: (I) The contact angle point position coordinates were obtained by the open-source algorithm. (II) The mineral labels of CTms in the 3×3 range centered on each contact angle point were counted and determined the mode. (III) The contact angle of each point was classified as the mineral category represented by the mode number. It is worth noting that the contact angle of pyrite is not analyzed because its contact angle number is too small.

2.4.3. Residual Oil Saturation and Pore Occupancy

The residual oil saturation was obtained by calculating the ratio of oil phase voxels to pore voxels based on each segmented data. In this paper, pore occupancy refers to the part of the pores occupied by oils, calculated by the ratio of oil phase volume and pore volume. It can be used to quantitatively analyze the relationship between pore size and fluid occupancy. The process is as follows: (I) The pore network model was extracted by the maximum ball algorithm developed by Raeini et al. [76] (see literature for details of the algorithm). The radius and center coordinates of all spheres (representing pore space) were obtained. (II) The center coordinates of all spheres were converted to pixel positions divided by image resolution and denoted as matrix Mball. (III) A sphere in matrix Mball was taken. In order to facilitate calculation, the minimum external cube of the sphere was denoted as the calculation region, and the number of pore voxels (num1) and the number of oil phase voxels (num2) in the calculation region were counted, respectively. Then, the pore occupancy was calculated by the ratio of num2 and num1 (shown in Figure 6). (IV) Step (III) was repeated for all spheres in the matrix Mball until the pore occupancy of all pores was calculated. Fifthly, the change rates of pore occupancy are calculated by
(2)Ri=POCoilPOCiPOCoil×100%i=1PV,20PV,50PV,
where R represents the change rate; POC represents the pore occupancy; and oil,1PV,20PV,50PV at the lower corner denotes saturated oil, flooding 1 PV, flooding 20 PV, and flooding 50 PV stage, respectively.

2.5. Distribution Pattern Classification of Residual Oil Based on Machine Learning

In previous studies, the shape factor and Euler number had been the most commonly used to divide the residual oil pattern (Yang et al., 2019; Song et al., 2020; An et al., 2016; Guo et al., 2018). In this study, to consider more shape features, the characteristic such as sphericity, shape factor, extent, axis rate of long to short, axis rate of middle to short, solidity, and Euler number of each residual oil block were calculated based on the regionprops3 function in Matlab software. The description and calculation method of each feature are as follows:

  • (i)

    Sphericity: Sphericity refers to the ratio of the surface area of an equal volume sphere to the actual surface area for each residual oil block. It is a parameter that characterizes the morphology of the residual oil blocks, and the closer to the block of the ball, the closer its sphericity is to 1

  • (ii)

    Shape factor: The shape factor is calculated by multiplying the volume by 6√π divided by the surface area to the power of 1.5, which can characterize the regularity degree of the residual oil blocks. The smaller the shape factor, the more irregular the blocks

  • (iii)

    Extent: The extent is the ratio of the volume of a residual oil block to the volume of the smallest cuboid containing the block, which can characterize the shape feature of residual oil blocks

  • (iv)

    Axis rate of long to short: Firstly, the length of the major axes of the ellipsoid that have the same normalized second central moments as the residual oil block was computed. Then, the axis rate of long to short can be calculated by the length ratio of the longest axes to the shortest axes

  • (v)

    Axis rate of middle to short: The axis rate of middle to short is the ratio of the length of the middle axes to the shortest axes. The axis rate of long to short and the axis rate of middle to short can jointly describe the extension characteristics of the residual oil blocks in three spatial directions

  • (vi)

    Solidity: The solidity refers to the proportion of a residual oil block’s volume to the volume of the smallest convex polygon that can contain the block, which can also characterize the shape feature of residual oil blocks

  • (vii)

    Euler number: The Euler number is an indicator of the connectivity of the 3D complex structure of residual oil blocks. It can be calculated by the number of connected components minus the number of holes then plus the number of enclosed cavities

In addition, artificial intelligence (AI) has been widely used in the oil industry due to its outstanding performance, including development effect evaluation [7779], data segmentation [8082], and flow simulation [83, 84]. In this study, the SOM machine learning method [55, 85] was used to automatically classify the residual oil pattern. The SOM method is a kind of artificial neural network with unsupervised learning and consists of treating the samples as vectors in an n-dimensional space, defined by a set of variables. The n-dimensional space will later be projected on a self-organized map with previously defined dimensions and several neurons. After competitive and collaborative steps, each original sample is represented by a best-matched unity (BMU). These BMUs are projected on maps preserving the topological relationships of similarity between the samples, measured by the Euclidean distance, aggregating the nonlinear relationships of the set of variables. The SOM method has been widely used in image classification [86, 87], environmental protection [88, 89], energy exploitation [55], etc. The SOM has two layers of network structure: the input layer and the output layer (also called the competition layer). In this study, for each residual oil block, seven feature parameters mentioned above were calculated as the train data. There were in total 15030 residual oil blocks in all displacement stages. So, the input layer size was set to 15030×7, i.e. there are 15030 neurons. The output layer size was set to 1×5, which means that the output layer has a one-dimensional structure composed of 5 neurons. After the training process, the 15030 inputs were mapped to 5 neurons in the output layer. Then, for each input, the winning neuron is determined based on its Euclidean distance from all neurons as a divided category. Finally, the 15030 residual oil blocks were automatically clustered into five patterns, as illustrated in Figure 7. During the training process, the neighborhood function was set to bubble, the activation distance was set to Euclidean, and the sigma value was set to 1.5. The initial learning rate was 0.5, and the training iterations number was set to 30,000. The characteristics of all residual oil patterns are as follows:

  • (i)

    Droplet: The droplet pattern oil dispersedly distributes in pore space. It is separated from larger residual oil blocks and can be deposited anywhere in the pore space, with small volume and regular in spherical or ellipsoidal

  • (ii)

    Film: The film pattern oil is mainly attached to the mineral surfaces. It has a small volume in the form of strips or cakes

  • (iii)

    Cluster: The cluster pattern oil distributes in several adjacent pore throats with narrowly connected throats

  • (iv)

    Network: The network pattern oil has a large area of contiguous composition. Compared with the cluster pattern, network pattern oil has a larger volume and better connectivity and only exists in the saturated oil stage

  • (v)

    Isolated: The isolated pattern oil is mainly located in a single pore. It is usually formed due to snap-off and trapped in the middle of the pores under the action of the Jiamin effect. It is irregularly shaped, approximately spherical, or ellipsoidal in appearance. It has a larger volume than the droplet pattern oil

3.1. Mineral Component Characteristics

Based on the QEMSCAN results (indicated in Figures 8(a)–8(b)), sample 1 contained 62.5% quartz, 16.2% feldspar, 3.0% siderite, 1.8% dolomite, 0.8% pyrite, 13.2% clay, and 2.5% other minerals. Sample 2 contained 57.3% quartz, 22.0% feldspar, 5.6% siderite, 1.7% dolomite, 0.4% pyrite, 10.1% clay, and 3.0% other minerals (shown in Figure 8(c)). The QEMSCAN results were further used to better identify and quantify the major minerals during the segmentation of CT images. Finally, the solid was segmented into quartz, feldspar, dolomite, siderite, clay, and pyrite, as illustrated in Figure 9. The proportion of each mineral and its surface calculated by CT data is displayed in Figure 10, where the quartz has the highest relative content, followed by feldspar and clay. Clay played a dominant role in the contribution of the surface, followed by quartz and feldspar. Therefore, the contribution of each mineral’s contribution to the mineral surface is not proportional to its relative content. Especially, clay provided about 50% of the mineral surface with about 10% content, making clay minerals have more opportunities to contact fluids because clay minerals have a larger surface area and higher surface roughness [90]. In addition, the content of quartz and feldspar was relatively high and therefore also provided considerable mineral surfaces. Thus, it can be inferred that clay, quartz, and feldspar significantly influence the multiphase flow transport mechanisms.

3.2. Wettability Variation of Different Mineral Surface

The subvolume of 700×700×700 voxels was extracted from the segmentation data of oil, 1 PV, 20 PV, and 50 PV stage, respectively, and the in situ contact angle was measured. The results are indicated in Figures 11 and 12. Figure 11 shows the contact angle histogram under different displacement stages. As can be seen, the contact angle is more widely distributed. In the saturated oil stage, sample 1 had a peak contact angle range of 65-70 degrees, different from the usual understanding. Generally, the reservoir’s wettability in the original state should be oil wettability [70, 91]. This difference may be due to no aging treatment after saturated oil. Then, at the 1 PV and 20 PV stages, the contact angle distribution shifted to the right and gradually became intermediate-wet. Although displaced by the brine, the sample had a long time of contact with oil under high temperature and high pressure, which played a role in aging. Under the comprehensive effect, the wettability shifted to intermediate. At the 50 PV stage, the oil phase volume decreased significantly, and the mineral surface was mainly contacted with brine, so the wettability shifted to water-wet again. The change of contact angle of sample 2 was similar to sample 1. The peak contact angle range was 60-65 degrees in the saturated oil stage, which was more hydrophilic than sample 1. Then in the 20 PV stage, the contact angle distribution shifted to the right. At the stage of 50 PV, the peak value range of the contact angle returned to 60-65 degrees. But sample 2 always displayed water-wet during the whole experiment with little change in contact angle.

Figure 12 demonstrates the contact angle histogram of each mineral surface. It can be seen from Figure 12(a) (sample 1) that the contact angle of siderite and dolomite changed significantly. The content of contact angle greater than 100 degrees and less than 30 degrees increased. The wettability of siderite and dolomite surfaces varies in opposite directions during the flooding, which is similar to the analysis results from the wettability-altered Ketton carbonate sample [62], due to the differential contact between the surface and fluids. The contact angle increased owing to the gradual aging of the surface in contact with crude oil, while the surfaces in contact with brine maintained their hydrophilic character. The contact angle of feldspar, quartz, and clay minerals showed an overall variation, which was shifted from water-wet to intermediate-wet under the comprehensive influence. However, after the 20 PV stage, with a substantial increase in brine phase volume, a large number of mineral surfaces resumed contact with brine. Until the 50 PV stage, the main peaks of contact angle of feldspar and clay minerals backed to 70 degrees and 72 degrees, respectively, showing water-wet characteristics. But, the main peak of quartz contact angle changed to 80 degrees, which remains in intermediate-wet.

Figure 12(b) (sample 2) illustrates that the wettability of siderite and dolomite surfaces also vary in opposite directions. When flooding to 50 PV, there was no apparent main peak in the contact angle distribution. The brine phase was fingered seriously along the hydrophilic and large pore-throat, resulting in a short contact time between the oil phase and mineral surface. So, feldspar, quartz, and clay minerals always displayed water-wet with little variation in contact angle. In summary, the wettability varied during the experiment, and each mineral displayed a different variation process. On the one hand, wettability evolution is affected by the properties of minerals [92]. For example, dolomite and calcite usually show stronger lipophilic properties. In contrast, quartz and feldspar are usually water-wetted [55, 71]. On the other hand, wettability is also affected by the fluid type and the contact time. Such as contact with crude oil for a long time, the wettability will vary to lipophilic by an adsorbed organic layer on surfaces [61]. If in contact with brine for a long time, it will show hydrophilic characteristics. In this study, some siderite and dolomite surfaces show lipophilic characteristics. Still, the content of such surface is low (only 2%, Figure 10), which has little effect on the overall wettability of the sample. Minerals such as feldspar, quartz, and clay provide the primary surface contribution and control wettability variation. According to the evolutionary process of contact angle, we can infer that sample 1 shows the weakly water-wet to intermediate-wet and sample 2 shows the water-wet. Different wettability will result in various impacts on the transport mechanisms of multiphase flow and distribution characteristics of fluids [67, 93, 94].

3.3. Variation of Residual Oil Saturation and Pore Occupancy

Figure 13 displays the variation of residual oil saturation (Sor). The initial Sor of sample 1 was 67.9%. Then, the Sor decreased rapidly (to 24.0%) until 1 PV stage, during which the oil phase volume decreased substantially. Then, the Sor decreased slowly (to 12.8%) until the 50 PV stage, and only a small amount of oil phase had been recovered. The final recovery rate was 81.1%. The Sor variation was different between sample 1 and sample 2. The initial Sor of sample 2 was 71.3%. After flooding to 1 PV, the Sor decreased rapidly to 44.3% and continued to 18.8% after flooding to 20 PV. However, the Sor remained unchanged until the 50 PV stage. The final recovery was 79.2%, lower than sample 1, indicating higher recovery under weakly water-wet to intermediate-wet sample.

The pore-throat distribution was statistically analyzed, as indicated in Figure 14. The pore-throat size distribution of the two samples was similar. Still, the ratio of pore in sample 2 with a radius greater than 50 μm, and the ratio of the throat with a radius greater than 40 μm was slightly higher. In addition, the average coordination number of samples 1 and 2 were 2.34 and 2.20, respectively, and the average pore-throat ratio was 1.40 and 1.31, respectively. Therefore, the pore structure of the two samples is very similar, which has a close effect on multiphase flow transport mechanisms.

Figure 15 displays the pore occupancy and its evolution. In sample 1 (Figure 15(a)), the pores with a radius of 90-100 μm had the highest occupancy rate after saturated oil. Because this part of the pores had a large pore radius, leading to small oil phase flow resistance and sufficient filling. The high capillary controlled the lowest occupancy of pores with a radius of 0-10 μm, which is resistant to oil flow. The pore occupancy in the radius of 20-50 μm also had a high level. It can be observed from the line graph that the oil phase volume proportion in a pore with a radius of 0-10 μm had the highest change rate until flooding to 1 PV and the lowest variation was found in the pore with a radius of 70-80 μm. The oil phase in the rest of the pores was displaced more uniformly under the weakly water-wet to intermediate-wet (Figure 11). The displacement front of the brine was relatively uniform (Figure 16), and the residual oil in the major pores was effectively displaced, so the Sor was significantly reduced during this period (Figure 13). When flooding to 20 PV, the oil phase volume proportion in the pores with a radius of 0-30 μm and 70-100 μm decreased further. The pore occupancy in the pores with a radius of 30-70 μm was relatively stable. The recovery rate was mainly contributed by the large pores (radius of 70-100 μm) and tiny pores (radius of 0-30 μm). When flooding to 50 PV, the oil phase volume in the pores with a radius of 30-70 μm began to effectively decrease. At the same time, the wettability shifted back to hydrophilic (Figure 11), with the extent of the brine spread range. It is worth noting that some oil phases re-entered the pore space with a radius of 0-20 μm and 80-100 μm and were trapped in the pore space, resulting in the addition of pore occupancy. In short, the displacement front is relatively uniform during the whole water flooding process (Figure 16), with continuous reduction of oil phase volume. The variation of oil phase volume proportion in each pore is relatively uniform. After 1 PV stage, the oil phase in the large and small pores is first displaced, followed by the oil phase in the intermediate pore.

In sample 2 (Figure 15(b)), the evolution of pore occupancy was similar to in sample 1 after saturated oil. The pores with a radius of 90-100 μm had the highest occupancy rate. In contrast, the pores with a radius of 0-10 μm had the lowest occupancy rate. The pore occupancy in the radius of 30-60 μm also had a high level. It can be seen from the line graph that the oil phase volume proportion in the pore with a radius of 60-100 μm decreased dramatically until flooding to 1 PV, which illustrated that the brine fingered along with the large pore space (Figure 16). The oil phase was displaced in large quantities, significantly contributing to oil recovery. In contrast, the change rate of oil phase volume proportion in pores with a radius smaller than 60 μm was low. The brine was preferred to flow along with the large pores with low flow resistance under the water-wet (Figure 11), and the front edge showed obvious fingering. When flooding to 20 PV, pores with a radius of 0-60 μm began to be spread, contributing to the recovery. Still, the pore occupancy in pores with a radius greater than 60 μm was stable. The oil phase was reinjected in pores with a radius of 70-80 μm from other locations, increasing the occupancy. At this time, each pore size has reached a high displacement efficiency. When flooding to 50 PV, occupancy in the pores with a radius of 0-30 μm had a significant reduction with the capillary force action. But, the content of this part of the crude oil was low, with less significance for recovery. And there is no obvious change of pore occupancy in pores with a radius greater than 30 um. In short, water flooding 0-20 PV is the main oil production period. After the 20 PV stage, only small pores have a slight recovery contribution. In water flooding, the water phase preferentially flows along with the large pores, with an obvious fingering phenomenon. Then, brine phase gradually extends to the small pores, and finally, only the oil in the small pores is displaced (Figure 16).

3.4. Changes in Residual Oil Patterns

As the analysis above, wettability significantly affects the multiphase flow transport mechanisms. This section further revealed the changes in residual oil patterns during the experiment process. All the residual oil blocks were divided into five patterns by the SOM method mentioned in Section 2.5, including droplets, films, clusters, networks, and isolated. Figure 17 displays the content variation of each residual oil pattern.

In sample 1, the residual oil content of the networks was dominant after saturated oil, and the oil phase was contiguous in an extensive range. After being injected into the sample, the brine advanced preferentially along the water film attached to the mineral surface, and the thickness of the water film gradually increased. At the narrow throat, the water film on both ends of the wall merged, i.e., the snap-off effect. Under the influence of the snap-off, the residual oil of the network pattern in a large contiguous area was gradually dispersed, forming the cluster and isolated oil. With water flooding, some cluster oil will be further snapped off and form more dispersed isolated oil, decreasing the cluster oil content. In contrast, isolated oil content continues to enhance. At the same time, due to long-term contact with crude oil, wettability inversion occurred on some mineral surfaces, forming film oil. Part of the cluster and isolated oil was trapped in the pore space due to the capillary force, and it was difficult to displace them without changing the development method (An et al., 2016). Finally, the residual oil was mainly in the form of the cluster and isolated with a small amount of film oil.

Sample 2 has the same distribution of the residual oil pattern after saturated oil as sample 1. After brine injection, the network oil was dispersed to form the cluster and isolated due to the snap-off effect. When flooding to 20 PV, some cluster oil will be further snapped off and form more dispersed isolated oil. But after the 20 PV stage, the content of cluster and isolated oil no longer changed significantly. In addition, when flooding to the 50 PV, the sample wettability reshifted to water-wet again (as described in Section 3.2). So the adsorption capacity of some mineral surfaces to the residual oil has reduced, which resulted in a decrease in the film oil content. Finally, the patterns of residual oil were also mainly cluster and isolated, with a small amount of film oil.

In this study, we carried out a new experimental procedure to reveal the impacts of the mineral surfaces’ wettability on fluids distribution characteristics and residual oil formation mechanisms during water flooding. We proposed a new mineral surface analysis algorithm to extract the mineral surface and classify its mineral categories. Each mineral in situ contact angle was measured to reveal the wettability variation. Pore occupancy in each pore was calculated to quantify the distribution characteristics of fluids. The self-organizing map (SOM) method was first used to classify the residual oil patterns.

Based on our analyses, we found that the clay minerals make the main contribution to the mineral surface, which provides about 50% of the mineral surface with about 10% content. The wettability of siderite and dolomite surfaces varies in opposite directions during the flooding. In contrast, the contact angles of feldspar, quartz, and clay minerals show an overall deviation related to the properties of mineral components and the contact time between solid and fluid. However, the total surface area of feldspar, quartz, and clay minerals occupies 95% of one of all the minerals, so these three minerals control the wettability of rock. The water flooding front of the weakly water-wet to intermediate-wet sample is uniform, and the residual oil in each pore is effectively displaced, resulting in high final recovery. In the water-wet sample, the water phase prefers to enter the large pores and form fingers and then gradually extends to tiny pores, with lower recovery. Under snap-off and wettability in version control, the residual oil distribution pattern evolves gradually from network to cluster, isolated, and film.

Due to the limitations of the subject confidentiality agreement, the data that support the findings of this study should be reasonably requested from the corresponding author.

The authors declare that there is no conflict of interest regarding the publication of this paper.

This work was funded by the National Natural Science Foundation of China (Grant No. U19B2006), the Support Program of China Postdoctoral Innovative Talents (Grant No. BX2021373), the Fundamental Research Funds for the Central Universities (Grant No. 22CX06006A), Postdoctoral Application Research Project of Qingdao City (qdyy20210088), the Major Scientific and Technological Projects of CNPC (Grant No. ZD2019-183-006), and Jidong Oilfield Science and Technology Special Project (Grant No. JDYT-2019-JS-134). Special Thanks are due to M.S. Rongji Zheng and M.S. Yiting Sun of China University of Petroleum (East China) for processing the data and revising the manuscript.

Exclusive Licensee GeoScienceWorld. Distributed under a Creative Commons Attribution License (CC BY 4.0).