Study on the Effects of Wettability and Pressure in Shale Matrix Nanopores Imbibition during Shut-in Process by Molecular Dynamics Simulations

1. Introduction

Shale reservoirs exhibit high tightness, low permeability, heterogeneity, and extensive nanopore development [1]. Conventional development techniques typically suffer from inefficiency. However, the utilization of multi-stage horizontal well volume fracturing technology for production purposes represents an effective approach [2,3]. The creation of stimulated fracture networks through hydraulic fracturing enhances the oil and gas seepage capacity of shale reservoirs. The oil present in the fractures formed during hydraulic fracturing can be extracted by utilizing the original formation pressure. Additionally, water imbibition into the matrix pores can displace shale oil present in the nanopores, enabling its extraction to the surface and facilitating EOR [4]. After the fracturing, it is common practice to shut the well for a certain period to allow for the full manifestation of the imbibition effect.

Spontaneous imbibition refers to the process in which porous media imbibe a wetting fluid and displace the original non-wetting fluid [5,6]. In low-permeability tight oil reservoirs, capillary imbibition serves as a significant mechanism for EOR, with capillary force acting as the primary driving force for spontaneous imbibition oil recovery [7]. Numerous scholars have conducted extensive experimental studies on the mechanism of imbibition oil recovery. The findings reveal that spontaneous imbibition within the shale matrix is influenced by various factors, including reservoir fluid, reservoir rock, and physical parameters such as pore structure, pore size distribution, oil-water interfacial tension, solid-liquid interface interaction, emulsification, wettability, temperature, shale oil viscosity, initial water saturation, and permeability, among others [8,9,10,11,12,13]. However, the wettability of the rock wall plays a crucial role in determining the effectiveness of spontaneous imbibition oil recovery. A more hydrophilic core leads to a faster spontaneous imbibition speed and higher recovery. Therefore, surfactants are frequently employed to improve the wettability of oil-wet reservoirs with limited imbibition oil recovery effects [14,15,16]. Secondly, due to the high rock fracture pressure associated with hydraulic fracturing, the pressure within fractures after fracturing is typically higher than the pressure within the matrix pores. As a result, both spontaneous imbibition within shale matrix pores and forced imbibition under the influence of high shut-in pressure take place. Thus, the shut-in pressure also serves as a critical factor that impacts the imbibition for EOR [17]. While traditional experimental methods can assess the imbibition oil recovery effect under various influencing factors, they have limitations in evaluating the oil recovery effects solely by calculating parameter changes from a macro perspective. Experimental research methods targeting various influencing factors include traditional experimental techniques such as weighing and volume observation methods, as well as advanced technologies, such as Nuclear Magnetic Resonance (NMR) and Computed Tomography (CT) technology [18,19,20,21]. Such experimental methods fail to directly depict the evolving characteristics of the oil and water phases within the pores and fall short of elucidating the microscopic mechanism underlying the imbibition oil displacement process.

With the rapid advancement of computer technology, MD simulation methods have gained widespread utilization in the field of oil and gas seepage [22,23]. MD simulation starts from the atomic and molecular scale and employs classical Newtonian mechanics to describe the movement of atoms, thereby capturing the dynamic characteristics of the modeled system [24]. The MD method has reached a relatively mature stage, characterized by minimal artificial influences during the simulation process and high credibility of the obtained results. Consequently, it has emerged as a crucial tool for unraveling the intricate mechanisms governing nanoscale fluid migration. Wang et al. [25] employed MD simulations to investigate the impact of nanoparticles on spontaneous water uptake in highly confined channels. They proposed a competitive mechanism for the spontaneous imbibition of nanofluids in capillaries by integrating the dynamic process of spontaneous imbibition, the contact angle of water in the capillary tube, and the relationship between displacement and time. Similarly, Yang et al. [26] conducted MD simulations to study the imbibition of octane and water into graphite and quartz nanopores. Their findings demonstrated that oil can be imbibed into dense organic nanopores of rocks at a faster rate than predicted by the classical imbibition model (Handy model). In another study, Nabin Kumar Karna et al. [27] performed large-scale atomistic simulations to investigate the capillary water imbibition phenomenon in pore silica nanochannels with varying heights (4-18 nm). They proposed an expression for the evolution of the dynamic contact angle in nanochannels over time, which, when combined with the Bosanquite equation, effectively explained the initial capillary rise. Additionally, Sang et al. [28] examined the imbibition of hydrocarbons in nano-kerogen pores. However, the aforementioned MD simulation studies on imbibition mainly focused on the imbibition of the wetting phase. They did not consider the displacement of the non-wetting phase from the pores. The two-phase imbibition mechanism, involving both the wetting and non-wetting phases, is more complex compared to the sole uptake of the wetting phase. Consequently, subsequent researchers conducted MD studies on the two-phase imbibition of oil and water, exploring various aspects such as the impact of polar components of shale oil on spontaneous imbibition [29] and the microscopic mechanism of water injection huff and puff for EOR [30]. Although some progress has been made, there is still limited research available on imbibition oil displacement, and the microscopic mechanism of the imbibition process involving both oil and water phases remains an area of active exploration.

In this study, we employed MD simulation to investigate the microscopic mechanism of imbibition oil displacement in matrix shale pores with different hydrophilicity and under different shut-in pressures. Initially, we visualized and quantitatively analyzed the spontaneous imbibition process to elucidate the microscopic advancement mechanism of the imbibition front and transportation characteristics of shale oil in different hydrophilic systems. Next, we analyzed the fluid migration characteristics during the forced imbibition process in different hydrophilic systems. Subsequently, based on the simulation results and certain theoretical derivations, we evaluated the synergistic effect of pressure difference and wettability on imbibition oil displacement.

2. Results

2.1. Shale oil migration in different hydrophilic nanopores during imbibition.

The spontaneous imbibition systems consisting of quartz walls W_A and W_B, along with oil-water, are named spontaneous imbibition system I (SI I) and spontaneous imbibition system II (SI II) respectively. Based on the static test results, symbols I-II represent the strongly hydrophilic system and weakly hydrophilic system respectively. During the spontaneous imbibition process, the fluid pressure surrounding the shale matrix is in equilibrium with the original matrix pore pressure[31]. Therefore, the pressure exerted on both the left and right He plates is maintained at 40MPa.

Firstly, MD simulations were performed on two spontaneous imbibition systems (SI I-II) for a duration of 1 ns. To intuitively reflect the behavior of oil-water phases in different hydrophilic systems during spontaneous imbibition, we calculated and presented the dynamic contact angle and simulated snapshots of imbibition systems I-II in Figure 1. Given the asynchronous advancement of the imbibition front on the upper and lower pore walls, we determine the dynamic contact angle as the arithmetic mean of the arithmetic means of the contact angles θ_d on the upper and θ_u lower pore walls. Meanwhile, H-bond exists if the distance of O−O is less than 3.5 Å and simultaneously the angle of H−O···O is less than 30° [32]. As shown in Figure 1a, as the wall wettability changes from strongly hydrophilic to weakly hydrophilic, the advancement distance of the imbibition front decreases significantly. Initially, water molecules located at the entrance penetrate the pore in the form of a meniscus. The progressive movement of the meniscus compels water molecules to continuously occupy nanopore space, leading to the displacement of the oil phase. As shown in Figure 1b, the dynamic contact angle in strongly hydrophilic system I is smaller than that in weakly hydrophilic system II and the dynamic contact angles in systems I and II are both less than 90°, which is in agreement with the results obtained from the static contact angle tests (Figure S5).

To intuitively reflect the axial migration behavior of the oil phase, we calculated the velocity distribution of the oil phase along the Z-axis and the two-dimensional (2-D) density distribution of the oil phase corresponding to the simulated snapshots of the imbibition systems I and II at 0, 500, and 1000 ps as shown in Figure 1a. The experimental test results show that the viscosity of crude decreases with increasing temperature. According to test results provided by the Xinjiang Research Institute of Petroleum Exploration and Development, Jimsar’s low sweet spot shale oil exhibits a high viscosity of tens to hundreds of mPa.S. However, the corresponding experimental data shows that at the reservoir temperature (315K), the viscosity of Jimsar shale oil has low viscosity close to that of water [33,34]. Therefore it can be assumed that the viscosity of shale oil under the reservoir conditions is equal to that of water. Furthermore, assuming the fluid density and viscosity do not vary spatially, The steady-state velocity profile for an incompressible laminar fluid confined between two quartz walls is parabolic and described by the classical Poiseuille equation.

$$v=-\frac{\Delta P}{2\eta }\left(z^2-\frac{w^2}{4}\right)$$

(1)

where ΔP is the pressure gradient along the flow direction (The pressure gradient ΔP during spontaneous imbibition is provided by capillary pressure); w and $\eta$ are the pore width and fluid (shale oil and water) viscosity respectively. Continuous hydrodynamics always assumes the streaming velocity of fluid vanishes at the interface. However, Recent studies have confirmed that pure hydrocarbons, such as octane and decane, exhibit a slip flow behavior when flowing on hydrophilic walls [35,36]. The slip length L_s is defined as the extrapolation distance to the surface location where the fluid velocity is equal to zero.

$$L_S=\pm \frac{v\left(z_{\mathrm{surf}}\right)}{{\left(\frac{dv}{dz}\right)}_{z_{\mathrm{surf}}}}$$

(2)

where $v\left(z_{\mathrm{surf}}\right)$ is the velocity of the fluid at the position of the fluid-solid interface z_surf. ${\left(\frac{dv}{dz}\right)}_{z_{\mathrm{surf}}}$ is the velocity gradient of the fluid at the fluid-solid interface.

Therefore, taking into account the effect of interface slip, the modified velocity distribution for the H-P equation is as follows:

$$v=-\frac{\Delta P}{2\eta }\left(z^2-\frac{w^2}{4}-w{L}_s\right)$$

(3)

Therefore, the L_s can be obtained by fitting the velocity distribution equation (3). Before extracting the velocity distribution curve, the Gibbs dividing surface (GDS) method was employed to determine the fluid-solid interface. However, considering that the imbibition front advances asynchronously on the upper and lower walls, for the accuracy and validity of the results, the velocity data on both sides of the center of the pore were fitted separately. Then the slip length is equal to the arithmetic mean of the slip lengths on the upper and lower walls of the pore. In the Figure 2a-b, the position where z = 0 represents the center of the pore.

From the velocity distribution diagram in Figure 2a,b, it can be observed that the actual velocity matches well with the fitted velocity and the oil phase velocity near the pore wall is much smaller. The occurrence of the slip phenomenon in the oil phase during the migration process proves that the overall push of the injected water to the oil phase is the main reason for the displacement of adsorbed oil molecules and that the oil phase has a larger slip length in more hydrophilic pores. The reason for the slip is that the strong attraction between hydrocarbons and quartz walls leads to the inability of water molecules to effectively peel off the oil phase. During the spontaneous imbibition process, as shown in Figure 2c, the shale oil still existing in the pores of either a strongly hydrophilic system or a weakly hydrophilic system always maintains good layered adsorption characteristics during spontaneous imbibition. Although the methyl and hydroxyl groups of the pore wall in the weakly hydrophilic system are uniformly arranged, the different forces of the methyl and hydroxyl groups on the hydrocarbons lead to the bending of the hydrocarbons that should be adsorbed parallel to the wall. However, the adsorbed oil in the pore walls of the weakly hydrophilic system has a higher adsorption strength (Figure S1), and the bending of the adsorbed hydrocarbons does not affect the final simulation results. To further elucidate the flow behavior of shale oil during spontaneous imbibition, the mass density distribution of shale oil and the number density distribution of different hydrocarbon components after 3ns of structural optimization, and the molecular number changes of different hydrocarbon components over time in the pores during spontaneous imbibition were calculated in Figure 3. It can be found, in agreement with the 2-D density distributions, that a distinct aggregation of hydrocarbon components on both sides of the pore wall occurs (Figure 3a,d). The oil in the pore contains three distinct adsorption layers, each with a width of 5 Å. Figure 3a,d showed that the average bulk density of shale oil calculated by MD simulation from z = -11 Å to z = 11 Å was 0.76 and 0.758 g/cm³ respectively, while the density measured by the experiment was 0.81 g/cm³ [33]. This was because C₂₀ was applied in the MD model to uniformly represent the C₂₀₊ molecules of the actual shale oil components, resulting in a lower density calculated by the model than that measured by the experiment. The oil within 5 Å of the quartz wall surface is defined as the first adsorption layer. The second and third adsorption layers are defined analogy. As shown in Figure 3b,e, the heavier hydrocarbon components generally have higher adsorption density in the first adsorption. There is no obvious peak density of C₄ in the first adsorption layer in the weakly hydrophilic system II because the pore wall modified by methyl groups of the weakly hydrophilic system II has a stronger adsorption effect on hydrocarbons. Therefore, the pore walls of the weakly hydrophilic system have a stronger adsorption effect on heavy and medium components, leading to a further reduction in the amount of C₄ adsorbed near the pore wall and a free state almost throughout the pore channels. During the imbibition period, the number of different hydrocarbon components in the first adsorption layer gradually decreases, and the number of hydrocarbons in the pores with more hydrophilic walls decreases faster (Figure 3c,f).

However, the reduction in the number of different hydrocarbon components in the first adsorption layer is by slip to the outside of the pore or by desorption into a bulk phase fluid. To evaluate the transport characteristics of different hydrocarbon components of adsorbed oil in pores, migration trajectories of C₄, C₈, and C₂₀ in the first adsorption layer after structural optimization during spontaneous imbibition are labeled as shown in Figure 4. The migration trajectories indicate that the heavy mass components in the adsorbed layer after structural optimization are always adsorbed parallel to the wall during the imbibition process. However, the desorption of light and medium components on the pore wall occurs as they are transported forward.

However, the phenomenon of a single molecule is not enough to reflect the overall motion characteristics of the adsorption layer. Therefore, the mean-squared displacement (MSD) in the X- and Z-directions of different hydrocarbon components in the first adsorption layer after structural optimization during spontaneous imbibition as shown in Figure 5 were calculated [37].

$$\mathrm{MSD}\left(t\right)=⟨{\left|{\mathit{r}}_i\left(t\right)-{\mathit{r}}_i\left(0\right)\right|}^2⟩$$

(4)

where ${\mathit{r}}_i\left(t\right)$ is the position of molecule i at time t and ${\mathit{r}}_i\left(0\right)$ is the initial position of molecule i.

For the hydrophilic system I, the MSD in the X-direction of the hydrocarbons is much larger than that in the Z-direction (Figure 5a,b), indicating that the first adsorbed layer of shale oil prefers to slip on the wall surface compared to desorption. As the hydrophilicity of the pore wall decreases, the MSD in the X-direction of the first adsorption layer of shale oil decreases and is similar to the MSD in the z-direction (Figure 5c,d), which means that the slip effect along the wall surface weakens and adsorbed oil molecules are neither easily desorbed nor diffused on the wall surface. Compared to the strongly hydrophilic system, in the weakly hydrophilic system, more adsorption sites are occupied by medium and heavy components, leading to slightly enhanced diffusion of light components, but there is no significant change in the MSD in the z-direction of the medium and heavy components (Figure 5a,c). Both in the strongly hydrophilic and weakly hydrophilic systems, the heavier hydrocarbon component has a smaller MSD in the Z-direction and a larger MSD in the X-direction, which suggests that the heavier component is more inclined to slip along the wall, but its slip effect is diminished. Therefore, the reason why the slip effect on the upper wall of the imbibition system I is more significant is because the adsorption effect of C₂₀ on the lower wall is more significant than on the upper wall (Figure 3b). The slip and diffusion of adsorbed oil on the wall surface are not conducive to imbibition and oil displacement. Therefore, how to effectively strip the adsorbed oil (mainly heavy components) during imbibition is a key issue to be considered to improve oil recovery.

2.2. Microscopic advancement mechanism of imbibition front

The driving force for spontaneous advancement of the meniscus at the molecular level arises from non-bonded interactions between the wall and water molecules, which encompasses both van der Waals and Coulomb interactions. The more hydrophilic walls have higher water-wall interactions and more hydrogen bonds are formed (Figure S2a,b). However, the approximately consistent change characteristics of Figure S2a and S2b indicate that Coulomb interaction predominantly governs this process due to the highly polar hydrophilic hydroxyl groups present on the pore wall. These hydroxyl groups can form hydrogen bonds with polar water molecules, resulting in a strong Coulombic interaction between water molecules and the pore wall [38].

Furthermore, as shown in Figure 2a,b, the oil phase has a larger slip length in more hydrophilic pores, however, the slip of shale oil on the wall surface is significantly reduced in the weakly hydrophilic system. Therefore, the wettability limitation of the slip length indicates that the advancement of the imbibition front is closely related to the fluid-solid interface behavior. Previous studies have provided little insight into the microscopic advancement of the imbibition front during the spontaneous imbibition process. Therefore, based on the relevant research on water imbibition into a capillary [39] and visualization results from the MD simulation of present spontaneous imbibition, we have constructed a schematic diagram depicting the advancement process of the imbibition front during spontaneous imbibition, as shown in Figure 6c. The water in the pore is divided into three parts: “adsorbed water molecules” near the pore wall, “bulk water molecules” in the center of the pore, and “interface water molecules” at the oil-water interface. Shale oil can be divided into “adsorbed oil molecules” near the pore wall and the “bulk phase oil molecules” in the center of the pore.

The “bulk water molecules” are almost not affected by the pore wall, but the “adsorbed water molecules” are firmly adsorbed on the pore wall due to the strong Coulomb interaction. Moreover, due to the strong Coulomb interaction of the solid pore wall, the interfacial water molecules move forward while approaching the pore wall, eventually becoming “adsorbed water molecules”. Therefore, the “interfacial water molecules” are continuously spread on the pore wall to realize the continuous advancement of the imbibition front and finally the displacement of shale oil out of the pore. To verify the imbibition front advancement process in Figure 6c, two “interfacial water molecules” (M_iw1 and M_iw2) and one adsorbed oil molecule (C₂₀H₄₂) (M_abo) in imbibition system I at the initial moment were marked (Figure 6a), and then their movement trajectory from 0 to 1000 ps was observed. The visualization trajectory shows that M_iw1 and M_iw2 were adsorbed on the pore wall after advancing a certain distance in the pore, eventually becoming immobile adsorbed water molecules. Secondly, the observation of the trajectory of M_abo reveals that it moves approximately horizontally to the right along the pore wall, which is consistent with the trajectory of C₂₀H₄₂ in Figure 4. The results show that the water molecules of the weakly hydrophilic system take a longer time to adsorb to the same wall position as M_iw2, which is the reason why the diffusion ability of shale oil of the weakly hydrophilic system is not significant.

To quantify the migration process of the three labeled molecules in Figure 6a, the X and Z coordinate values of M_iw1, M_iw2, and M_abo at different moments were counted, as shown in Figure 6d,e. For M_iw1, although its X and Z coordinate values fluctuate locally, its X coordinate in general appears first gradually to increase and then stabilize, and its Z coordinate generally appears to gradually decrease and then stabilize, indicating that M_iw1 flowed along the imbibition direction (X direction) while approaching and adsorbed on the pore wall. The trajectory of M_iw2 is similar to that of M_iw1, but the X and Z coordinate values are stabilized more quickly, suggesting that it adsorbs to the pore wall more quickly compared to M_iw1. The Z coordinate value of M_abo fluctuates near the pore wall, and the X coordinate value continues to increase overall, which proves that M_abo is continuously displaced along the pore wall to the outside of the pore. The coordinate values of M_iw1, M_iw2, and M_abo are consistent with the trajectory in Figure 6a. Thus, the Microscopic advancement mechanism of the imbibition front during spontaneous imbibition is the competitive adsorption between “interfacial water molecules” and “adsorbed oil molecules”. To determine the reason for the slower diffusion of shale oil at the pore wall in the weakly hydrophilic system, an “interfacial water molecule” in the weakly hydrophilic system was labeled in Figure 6b, which has the same initial position and the same adsorption site on the wall as M_iw2 during the imbibition process. A comparison of Figure 6a and b reveals that it takes longer for “interfacial water molecules” in a weakly hydrophilic system to adsorb to the same position as M_iw2, which is the reason why adsorbed oils in the weakly hydrophilic system do not have a significant ability to diffuse on the wall.

Because of the asynchronous advancement of the imbibition front, the spontaneous imbibition process in the nanopore can be divided into two stages based on changes in replacement efficiency[40]:

$$displacementefficiency=\frac{N_{Os}-N_{Ns}}{N_{Os}}$$

(5)

where N_Os represent the number of oil molecules originally in the pore; N_Ns represent the number of oil molecules in the pore at time t.

As shown in Figure 7d, we observe the displacement efficiency increases approximately linearly in the first stage I (At 0-1621 ps). However, the increase in displacement efficiency in stages II (1621-2519 ps) slows down because the imbibition front reaches the end on the upper walls at 1621ps so that the competitive adsorption of "interfacial water molecules" of the imbibition front and "adsorbed oil molecules" on the pore wall only occurs on the lower wall. The displacement efficiency in stage I amounts to 69.23%, with a displacement efficiency per unit time of 0.043%ps^-1. Moreover, the displacement efficiency in stage II amounts to 30.77%, with a displacement efficiency per unit time of 0.034%ps^-1. That is to say, the displacement velocity in the second stage dropped to 0.77 times that in the first stage. The change in displacement efficiency also proves that spontaneous imbibition is the result of the competitive adsorption between "interfacial water molecules" and "adsorbed oil molecules".

In addition, the velocity distribution of water molecules in the X-direction in both spontaneous imbibition systems I and II at 500ps in Figure 7a,b show significant fluctuations at the oil-water interface. However, the velocity fluctuations of water molecules in the rest of the pore were relatively small. Therefore, the transport of bulk-phase water molecules in the pore is less affected by the wall surface and can be approximated as horizontal advancement. In other words, the essence of spontaneous imbibition oil displacement is also that water molecules outside the pore continuously adsorb and accumulate on the pore wall surface (Figure 7c).

Therefore, the essence of spontaneous oil displacement is the adsorption and accumulation of water molecules outside the matrix pores around the hydroxyl groups on the pore wall, which can be described by the radial distribution function g(r). Figure 8a,b depicts the g(r) between water molecules and hydroxyl groups on the pore wall. The first peaks of the water molecules in both systems I and II were about 0.275 nm away from the hydroxyl groups on the wall, which was close to the distance of the hydrogen bonds [41]. It proves that there are strong hydrogen bonding interactions between the water molecules and the hydroxyl groups on the wall and the water molecules gather and adsorb around hydroxyl groups because of the hydrogen bond interactions. The peak of g(r) and the area under the curve increase with time, suggesting a continuous rise in the coordination number of water molecules surrounding the hydroxyl groups on the pore walls, ie., water molecules spread forward along the pore walls. By comparing the g(r) values in Figure 8a and b, it can be found that the adsorption and aggregation of water molecules on the wall in the weakly hydrophilic system are weakened because of the decrease in the density of hydroxyl groups on the wall. The g(r) between oil and water, displayed in Figure 8c-d, is found to be small and remains stable over time, indicating the oil-water interface is stable and there is no significant miscibility between them. Thus, the displacement of oil occurs in a piston-like way. A clear oil-water interface plays a crucial role in preserving optimal interfacial tension, thereby ensuring an effective capillary pressure that serves as the driving force for spontaneous imbibition.

2.3. Forced imbibition under pressure difference

Spontaneous imbibition generally occurs under forced pressure (the difference between shut-in pressure and original pore pressure) during the shut-in period [42]. To simulate the forced imbibition process when the pressure difference is 5, 10, and 20 MPa, we apply a matrix pore pressure of 40 MPa to the He plate 2, and different shut-in pressures of 45, 50, and 60 MPa to the He plate 1. The forced imbibition systems corresponding to systems I and II are named FI I and FI II respectively.

The influence of shut-in pressure on imbibition during the shut-in period is vividly described through a series of snapshots at different moments (0, 250, 500, 750, 1000 ps) under varying pressure differences as illustrated in Figure 9a-b. These snapshots are compared with the corresponding spontaneous imbibition simulation snapshots. It can be intuitively observed that as the shut-in pressure increases, the imbibition effect also increases accordingly. Notably, the 2-D density distribution shows that the adsorbed oil still present in the pores remains consistently better-adsorbed structures during forced imbibition under pressure difference (Figure 9c).

To quantitatively characterize the stability of the adsorption layer during forced imbibition under pressure difference, the MSD in the Z-direction of different hydrocarbon components in the first adsorption layer after structural optimization during forced imbibition was calculated, as shown in Figure 10. In forced imbibition under pressure difference, the diffusion performance of the light component changes relatively more, and the diffusion ability of the medium component in the Z-direction changes relatively less (Figure 10a,b,d,e). The MSD of the heavy components in the Z-direction under the effect of pressure difference is not significant (Figure 10c,f), which means that the pressure difference has no significant desorption effect on the heavy components. Therefore, the axial pressure difference leads to changes irregularly in adsorption diffusion properties but does not effectively strip the adsorbed layer of shale oil (especially the heavy component).

Further observation of the imbibition front of the simulated snapshots of forced imbibition in Figure 9a-b shows that with the increase of pressure difference, the degree of curvature of the meniscus at the imbibition front tends to weaken (ie., The value of θ_d increases with the increase of shut-in pressure), and it is especially obvious at 20 MPa. In weakly hydrophilic system II, the phenomenon of boundary water lagging behind the water in the bulk phase was observed even under the pressure difference of 20 MPa, which means that the morphology of the meniscus at the imbibition front changes from concave to convex on the aqueous side. This may be due to wetting hysteresis caused by differential pressure [43].

Because the uneven distribution of mixed crude oil results in asynchronous advancement of the imbibition front and difficulty in accurately measuring the dynamic contact angle, other methods are needed to quantitatively characterize the wetting hysteresis phenomenon. The wetting hysteresis phenomenon is closely related to the velocity distribution of the fluid of the nanopore in the Z-direction. Thus, the slip velocity and peak velocity variations can quantitatively characterize the wetting hysteresis phenomenon and can be obtained by fitting equation (3), As shown in Figure 11a. Considering the asynchronous advancement of the imbibition front and the accuracy of the results, both the slip and peak velocities are obtained by arithmetic averaging, as shown in equations (6) and(7).

$$v_{slip}=\frac{v_{slip}^{lower}+v_{slip}^{upper}}{2}$$

(6)

$$v_{peak}=\frac{v_{peak}^{lower}+v_{peak}^{upper}}{2}$$

(7)

where $v_{slip}^{lower}$, $v_{slip}^{upper}$, and $v_{slip}$ are the slip velocities of the shale oil on the lower and upper walls, and their arithmetic mean respectively; $v_{peak}^{lower}$, $v_{peak}^{upper}$, and $v_{peak}$ are the peak velocities of shale oil on the lower and uppper sides of the center of the pore, and their arithmetic mean respectively. The fitting results for the systems I and II are presented in Figure 11b and 11c respectively.

Linear fitting again of the velocity obtained by fitting Equation (3) can obtain the slopes K_peak and K_slip respectively representing the speed changes of peak velocity and slip velocity with pressure difference. Both peak and slip velocities increase faster for the strongly hydrophilic system I than for the weakly hydrophilic system II, representing a more significant promotion of the imbibition in the strongly hydrophilic system by the pressure difference. In both strongly and weakly hydrophilic systems, K_peak is greater than K_slip, leading to a decrease in the dynamic contact angle (i.e., wetting hysteresis). Further comparison of the ratios of K_peak and K_slip (i.e., $\frac{K_{peak}}{Kslip}$) for the different imbibition systems, respectively, the ratio $\frac{K_{peak}}{Kslip}$ of the strongly hydrophilic system I is greater than that of the weakly hydrophilic system II, indicating that the difference between the peak and slip velocity in the strongly hydrophilic system changes less than that of the weak hydrophilic system, which means that the wetting hysteresis phenomenon of the weak hydrophilic system is more significant.

To further clarify the wettability and pressure difference correlation because of the wetting hysteresis, We have carried out certain theoretical derivations and analyses. Based on the reference [44], when forced imbibition occurs after shut-in, the equilibrium of forces on the fluid in the confined capillary can be expressed as:

$$F_c+F_p+f+F_i+F_g=0$$

(8)

where F_c is the capillary force; ΔF is the displacement force. f is the viscous force generated by the friction between fluid molecules and between fluid molecules and the solid pore wall; F_i is the inertial force; $F_g$ is gravity; Because the matrix pore width of tight shale reservoirs is extremely small and much smaller than the imbibition height, the influence of $F_g$ can be ignored [45]. Secondly, for pores with extremely small radii, viscous force $f$ dominates and the inertial force $F_i$ can be ignored [46]. Then equation (8) can be simplified to:

$$F_c+F_p+f=0$$

(9)

Assuming that the width between two walls, length in the Y-direction, and length in the X-direction of the slitpore are w, l, and L respectively, It is possible to convert $F_p$ into wl∆P. According to the Young-Laplace equation [47], the F_c in a single pore is:

$$F_c=\frac{2\sigma cos\theta }{w}wl=2\sigma l\cos \theta$$

(10)

where θ is the contact angle of the oil-water-wall three-phase system, which is assumed to be equal to the static contact angle and is not affected by other factors; σ is the oil-water interfacial tension. The total f of oil and water in a single pore is:

$$f=2xl{\tau }_w+2l\left(L-x\right){\tau }_o$$

(11)

where x is the position of the oil-water interface in the X-direction; τ_w and τ_o are the shear stress of water and oil respectively, which can be expressed as [48]:

$$\begin{array}{c}{\tau }_w=-\frac{12{\mu }_w}{w}v\\ {\tau }_o=-\frac{12{\mu }_o}{w}v\end{array}$$

(12)

where v is the imbibition speed, μ_w and μ_o are the viscosity of water and oil respectively. After bringing equation (10-12) into equation (9) and rewriting v as the derivative dx/dt of x with respect to time t, equation (9) is transformed into:

$$∆P{w}^2+2\mathit{\sigma w}\cos \theta -24\left[x{\mu }_w+\left(L-x\right){\mu }_o\right]\frac{dx}{dt}=0$$

(13)

After integrating equation (13) under the initial condition of x = 0 and t = 0, we get:

$$12\left({\mu }_w-{\mu }_o\right)x^2+24{\mu }_{0}Lx-\left(\mathsf{\Delta }\mathrm{P}{w}^2+2\sigma wcos\theta \right)t=0$$

(14)

Additionally, it should be noted that the shale oil model used in this study does not account for components such as C₂₀₊ and asphaltenes, implying that its viscosity could be lower under reservoir conditions. As previously mentioned, the oil displacement efficiency primarily exhibits a linear change during the imbibition process before the imbibition front reaches the outlet. Considering the low viscosity of the Jimsar shale oil under reservoir conditions, for simplification, the viscosity of shale oil under reservoir conditions can be approximated as μ = μ_w = μ_o = 0.36 mPa·s, ie., the viscosity of water under the same conditions obtained from the NIST database. Consequently, equation (14) can be simplified as follows:

$$24\mu Lx-\left(\Delta p{w}^2+2\sigma wcos\theta \right)t=0$$

(15)

Then the imbibition velocity at pressure difference is

$$v=\frac{x}{t}=\frac{\left(\mathsf{\Delta }P{w}^2+2\sigma w\mathrm{c}\mathrm{o}\mathrm{s}\theta \right)}{24\mu L}$$

(16)

The pressure differences are equal to 0, 5, 10, and 20MPa respectively in this study. When the pressure difference is equal to 0, Equation (16) is converted to v = σwcosθ⁄12uL, indicating that the spontaneous imbibition velocity is constant during spontaneous imbibition under reservoir conditions. This may be the reason why the oil displacement efficiency changes approximately linearly under reservoir conditions before the imbibition front reaches the outlet.

The imbibition velocity when the pressure difference is equal to 0 is called spontaneous imbibition velocity (v_SI), and the imbibition velocity when the pressure difference is equal to 5, 10, and 20MPa is called forced imbibition velocity. Assuming that when capillary imbibition occurs, the dynamic contact angle does not change with pressure difference, the ratio (R) of forced imbibition to spontaneous imbibition velocity at pressure difference is as follows:

$$R=\frac{v_{FI}}{v_{SI}}=\frac{\Delta P{w}^2}{2\sigma wcos\theta }+1$$

(17)

From equation (17), assuming that the three-phase contact angle $\theta$ is not affected by the pressure difference and is always equal to the static contact angle, the R of the weakly hydrophilic system should be greater than that of the strongly hydrophilic system at different pressure differences. However, the R obtained through MD simulations as shown in Table 1 shows that the R of the strongly hydrophilic system is larger than that of the weakly hydrophilic system, which further proves that the pressure leads to the occurrence of the wetting hysteresis phenomenon.

Some parameters of Equation (16) are: The interfacial tension σ under reservoir conditions obtained by MD simulation is 55.44mN/m (Figure S3) (The results indicate a good match between the simulation and experimental value of 53.2 mN/m); The static contact angles θ for the strongly and weakly hydrophilic walls are 53.9° and 78.97° respectively; The length L and width w of the pore is 15 and 6 nm. The imbibition speed at different pressure differences obtained by substituting the value of these parameters mentioned above into Equation (16) is renamed as the theoretical imbibition speed $v_{th}$ under different pressure differences. Meanwhile, we rename the imbibition velocity obtained through MD simulation in Table 1 as $v_{MD}$. The calculation method of the imbibition velocity obtained through MD simulation is: taking the average value of the fitted velocity curve in Figure 11a. Although the theoretically calculated and MD simulated values do not necessarily correspond to the real imbition velocity, their change characteristics can quantitatively reflect a certain imbibition mechanism.

Figure 12a,b illustrates the values of $v_{MD}$ and $v_{th}$, as well as the difference (Δv) (∆v=$v_{MD}$- $v_{th}$) between the them at different pressure differences. The results show basically that $v_{MD}$ and $v_{th}$ exhibit the same order of magnitude, confirming the reliability of the simulation outcomes. In imbibition systems I and II, as the pressure difference increases, both $v_{MD}$ and $v_{th}$ gradually increase, while Δv continues to decrease. The results demonstrate that an increase in pressure difference can enhance the imbibition effect, while also leading to a reduction in the effective imbibition velocity relative to the theoretical value because of the wetting hysteresis phenomenon. Moreover, ∆v decreases more rapidly in the weakly hydrophilic system relative to the strongly hydrophilic system, suggesting that a more pronounced wetting hysteresis phenomenon leads to a more pronounced loss of the imbibition velocity relative to theoretical values. The loss of imbibition velocity relative to the theoretical velocity is caused by the reduction of capillary pressure due to wetting hysteresis. Consequently, excessively high imbibition pressure impedes the maximization of economic benefits and reduces the underutilization of natural power.

3. Models and Methodology

3.1. Model system

All molecular models were constructed through Materials Studio software developed by Accelrys Company in the United States. The molecular model constructed in this article mainly includes solid quartz walls, nanopores composed of quartz walls, and fluid components (shale oil and water).

Scholars [49,50] have conducted mineral components analysis on shale core samples from the Jimsar reservoir by using a device of X-ray diffraction. The results show that the mineral composition of the Lucaogou Formation shale matrix is complex and diverse, but quartz has the highest relative content and plays an important role. Quartz is often used to simulate inorganic nanopores[40,51,52,53]. Therefore, quartz is used in this work to represent shale reservoir matrix pores. The initial quartz unit cell SiO₂ Quartz was taken from the structure database of Materials Studio software and then was cleaved along the ($0\overline{01}$) crystal face. After being cleaved, the crystal face was replicated to get the final quartz wall. The complexity and randomness of the mineral distribution in shale reservoirs lead to a heterogeneous distribution of wettability as well [54]. To characterize the uneven wettability of shale reservoirs, silica was modified concerning existing research methods. The hydrophilic modification of the wall surface (named W_A) was first realized by adding hydrophilic group hydroxyl (-OH) to all chemically unsaturated Si atoms of the ($0\overline{01}$) crystal surface (surface in contact with oil), with a surface hydroxyl density of 9.6 nm⁻² which was consistent with the results of the crystal chemistry calculations (5.9–18.8 nm⁻² ) [55]. Then hydroxyl hydrogens were substituted with a certain amount of hydrophobic group methyl (-CH₃) to obtain the other wall surface (named W_B) with different wettability. To make the wall surface as homogeneous as possible, all the hydroxyl and methyl groups were evenly distributed on the wall surface of the quartz. The hydroxyl and methyl density on the surface in contact with the oil of the modified quartz wall was 7.2 nm⁻² and 2.4 nm^-2. The back surface and left and right end surfaces of the quartz walls were all fully hydroxylated. Considering that spontaneous imbibition in oil-wet pores cannot occur effectively, only hydrophilic matrix pores of the reservoir are considered in this paper. Each wall has a thickness of 21.6 Å and dimensions in the X and Y directions of 150Å and 29.5 Å, respectively. Quartz walls with different hydrophilicity are presented in Figure S4.

Before conducting MD imbibition simulations, the static wettability of different walls was first quantitatively characterized (The corresponding results are presented in Figure S5). Then on this basis, the molecular model of imbibition oil displacement was constructed as shown in Figure 1. The nanopore of the imbibition model consisted of two quartz walls with the same modification as shown in Figure S4. The two walls of the pore were placed parallel to the X-Y plane, and the vertical distance between their wall surfaces in the Z direction was 60 Å. For the shale oil in the pore, considering the complexity of shale oil components, a single oil component cannot truly reflect the occurrence and flow state of shale oil in the pore. According to Table S2, in this study, a ternary system composed of n-Butane (C₄H₁₀), n-Octane (n-C₈H₁₈), n-Eicosane (n-C₂₀H₄₂) was selected to represent shale oil in the lower “sweet spot” of the Lucaogou Formation in the Jimsar Sag, Xinjiang oilfield. n-Butane (n-C₄H₁₀), n-Octane (n-C₈H₁₈), and n-Eicosane (n-C₂₀H₄₂) respectively represent the light, medium and heavy components in Jimusar shale oil. The molar fractions and number of molecules of each hydrocarbon component of the shale oil model in the pore are shown in Table S2. Water molecules representing the fracture, i.e., fluid surrounding the matrix, are placed to the left of the pore. Before the imbibition oil displacement simulation, the shale oil was fully adsorbed on the pore wall through 3ns equilibrium molecular dynamics (EMD) at reservoir conditions (353.15K and 40MPa), thereby obtaining a reasonable shale oil density. During the imbibition simulation, He plates were placed on the left and right sides of the oil-water-wall system, and the imbibition oil displacement process was controlled by controlling the pressure attached to the He plate. In this work, it is assumed that the shut-in pressure is equal to the fluid pressure around the matrix. The pressure on He plate 2 always maintains the matrix pore pressure (P_m) of 40MPa, while the pressure exerted on He plate 1 represents the shut-in pressure (P_s). When the pressure on He plate 1 is equal to a matrix pore pressure of 40MPa, water is spontaneously imbibed into the pore. When the pressure applied to He plate 1 exceeds 40MPa, it is the process of forced imbibition oil displacement. The simulation time of imbibition oil displacement is 1ns.

3.2. Simulation Details

Different force fields have different characteristics and applicable conditions. Therefore, choosing an appropriate force field for a specific research system is key to the accuracy and reliability of MD simulation results. In this paper, we use different force field potential energy parameters for different molecular groups. The ClayFF force field was developed by Randall T et al [56] for hydrates and multi-component minerals and their interaction with fluid interfaces and has been successfully applied to hydroxylated silica surface-fluid systems [57,58,59,60]. This paper uses the ClayFF force field to describe the silica crystal and surface hydroxyl groups. Referring to the research of Aleksandr Abramov, the remaining surface methyl groups are described by the DREIDING force field [61]. The OPLS-AA force field [62] is used to describe n-hydrocarbons in shale oil. The OPLS force field can accurately describe hydrocarbon molecules’ thermodynamic properties and structural characteristics [63]. The SPC/E potential energy model [64] is used to describe water molecules. SPC/E can reproduce the dynamic physical behavior of water molecules on the interface and has good compatibility with the CLAYF force field and OPLS-AA force field [65,66,67,68]. The helium plate and quartz base plate are considered rigid bodies. Periodic boundary conditions are applied in all directions. This article uses Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) open-source program software to perform MD simulations [69], with a time step of 1fs. The simulation results were visualized using VMD software [70]. Intermolecular non-bonded interactions (Enonbonded) can be described by the standard 12-6 Lennard-Jones (LJ) potential (E_LJ ) [54] and Coulomb potential respectively (E_coulomb) [58,71].

$$E_{nonbonded}=E_{LJ}+E_{coulomb}=4{\epsilon }_{ij}\left[{\left(\frac{{\sigma }_{ij}}{r_{\mathrm{i}\mathrm{j}}}\right)}^{12}-{\left(\frac{{\sigma }_{ij}}{r_{\mathrm{i}\mathrm{j}}}\right)}^6\right]+\frac{q_i{q}_j}{{4\pi \epsilon }_0{r}_{ij}}$$

(18)

where r_ij stands for the distance between two atoms i and j; εij and σij represent the LJ well depth and the zero-potential distance between atoms i and j respectively; ε₀ represent the dielectric constant of a vacuum; And q_i and q_j represent the charge of atoms i and j respectively. The interactions between different atoms are calculated using the Lorentz-Berthelot mixing criterion:

This study uses MD simulation to elucidate the microscopic mechanism of water imbibition and oil displacement in the different wettability pores of the shale matrix under different shut-in pressures. The main conclusions are as follows:

$${\sigma }_{ij}=\frac{1}{2}\left({\sigma }_{ii}+{\sigma }_{jj}\right)$$

(19)

$${\epsilon }_{ij}=\sqrt{{\epsilon }_{ii}{\epsilon }_{jj}}$$

(20)

4. Conclusions

This study uses MD simulation to elucidate the microscopic mechanism of water imbibition and oil displacement in the different hydrophilic pores of the shale matrix under different shut-in pressures. The main conclusions are as follows:

(1) The slip phenomenon suggests that the strong attraction between hydrocarbons and quartz walls leads to the inability of water molecules to effectively peel off the oil phase. The heavy hydrocarbon component in the adsorbed layer tends to slip on the walls and the light component tends to desorb from the pore wall. However, the heavier hydrocarbon components have poorer slip ability on the wall. Weak wall hydrophilicity and slip of shale oil on the wall surface can significantly reduce the efficiency of imbibition displacement. How to effectively strip the adsorbed oil (especially heavy components) during imbibition is a key issue to be considered to improve oil recovery.

(2) The displacement of oil occurs in a piston-like way during spontaneous imbibition. The microscopic advancement mechanism of the imbibition front is the competitive adsorption between “interfacial water molecules” at the imbibition front and “adsorbed oil molecules” on the pore wall. The essence of spontaneous imbibition is the adsorption and aggregation of water molecules onto the hydroxyl groups present on the pore wall. The better hydrophilicity of the matrix pore's wall facilitates higher levels of adsorption and accumulation of water molecules on the pore wall and requires less time for “interfacial water molecules” to compete with “adsorbed oil molecules” because interfacial water molecules have a stronger competitive adsorption effect on the more hydrophilic pore.

(3) During the forced imbibition process, as the pressure difference increases, the imbibition effect gradually increases. The pressure difference acts mainly on the bulk oil thus causing wetting hysteresis to occur. Meanwhile, shale oil still existing in the pore always maintains a good stratified adsorption structure. During forced imbibition, the actual capillary pressure gradually decreases because of the wetting hysteresis phenomenon and there is a loss in the imbibition velocity relative to the theoretical value. Simultaneously, the decline in hydrophilicity further weakens the synergistic effect of pressure difference because of the more pronounced wetting hysteresis. Consequently, excessively high shut-in pressure impedes the maximization of economic benefits and reduces the underutilization of natural power.

The study presented in this article investigates the microscopic mechanism and process of imbibition oil displacement in matrix shale pores under different wettabilities and shut-in well pressures. The findings offer valuable insights into enhancing imbibition oil recovery methods, optimizing well stew time, and improving production in shale reservoirs.