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

Shut-in after fracturing is generally adopted for wells in shale oil reservoirs, and imbibition occurring in matrix nanopores has been proven as an effective way to improve recovery. In this research, a molecular dynamics (MD) simulation was used to investigate the effects of wettability and pressure on nanopore imbibition during shut-in for a typical shale reservoir, Jimsar. The results indicate that 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 involves the adsorption and aggregation of water molecules onto the hydroxyl groups on the pore wall. The flow characteristics of shale oil suggest that the overall push of the injected water to the oil phase is the main reason for the displacement of adsorbed oil molecules. Thus, shale oil, especially the heavy hydrocarbon component in the adsorbed layer, tends to slip on the walls. However, the weak slip ability of heavy components on the wall surface is an important reason that restricts the displacement efficiency of shale oil during spontaneous imbibition. The effectiveness of spontaneous imbibition is strongly dependent on the hydrophilicity of the matrix pore’s wall. The better hydrophilicity of the matrix pore 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. During the forced imbibition process, the pressure difference acts on both the bulk oil and the boundary adsorption oil, but mainly on the bulk oil, which leads to the occurrence of wetting hysteresis. Meanwhile, shale oil still existing in the pore always maintains a good, stratified adsorption structure. Because of the wetting hysteresis phenomenon, as the pressure difference increases, the imbibition effect gradually increases, but the actual capillary pressure gradually decreases and there is a loss in the imbibition velocity relative to the theoretical value. Simultaneously, the decline in hydrophilicity further weakens the synergistic effect on the imbibition of the pressure difference because of the more pronounced wetting hysteresis. Thus, selecting an appropriate well pressure enables cost savings and maximizes the utilization of the formation’s natural power for enhanced oil recovery (EOR).


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 the 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 velocity 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 the 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 during shut-in for the shale reservoir.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.

Shale Oil Migration in Different Hydrophilic Nanopores during Spontaneous Imbibition
The spontaneous imbibition systems consisting of quartz walls W A and W B , along with oil-water, are named the 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 40 MPa.
Firstly, MD simulations were performed on two spontaneous imbibition systems (SI I and SI 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 (θ d ) and simulated snapshots of SI I and SI II in Figure 1.Given the asynchronous advancement of the imbibition front on the upper and lower pore walls, we determine the θ d as the arithmetic mean of the contact angles (θ u ) on the upper pore wall and (θ l ) on the lower pore wall (the inset in Figure 1b).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 the nanopore space, leading to the displacement of the oil phase.As shown in Figure 1b, the dynamic contact angle in the strongly hydrophilic system is smaller than that in the weakly hydrophilic system, and the θ d in SI I and SI II are both less than 90 • , which is in agreement with the results obtained from the static contact angle tests (Figure S5d).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 SI I and SI II at 0, 500, and 1000 ps, as shown in Figure 2. The experimental test results show that the viscosity of crude oil 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 (315 K), the viscosity of Jimsar shale oil has a low viscosity close to that of water [32,33].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 Hagen-Poiseuille (H-P) equation [34].
where ∆ is the pressure gradient along the flow direction (the pressure gradient ∆ during spontaneous imbibition is provided by capillary pressure);  and  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 [34,35].The slip length ( ) is defined as the extrapolation distance to the surface location where the fluid velocity is equal to zero [36].
where ( surf ) is the velocity of the fluid at the position of the fluid-solid interface ( ).
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: 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 SI I and SI II at 0, 500, and 1000 ps, as shown in Figure 2. The experimental test results show that the viscosity of crude oil 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 (315 K), the viscosity of Jimsar shale oil has a low viscosity close to that of water [32,33].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 Hagen-Poiseuille (H-P) equation [34].
where ∆P is the pressure gradient along the flow direction (the pressure gradient ∆P during spontaneous imbibition is provided by capillary pressure); w and η 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 [34,35].The slip length (L S ) is defined as the extrapolation distance to the surface location where the fluid velocity is equal to zero [36].
where v(z surf ) is the velocity of the fluid at the position of the fluid-solid interface (z surf ).Therefore, taking into account the effect of interface slip, the modified velocity distribution for the H-P equation is as follows: Therefore, the L S can be obtained by fitting the velocity distribution of shale oil in Figure 2a,b using Equation (3).Before extracting the velocity distribution curve, the Gibbs dividing surface (GDS) method was employed to determine the fluid-solid interface [37].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 L S is equal to the arithmetic mean of the L S on the upper and lower walls of the pore.As shown in the inset in Figure 2a, 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 velocity distribution matches well with Equation (3), 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 L S 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.As shown in Figure 2c, the shale oil still existing in the pores of both a strongly hydrophilic system and 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 3 ns 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 about 5 Å. Figure 3a,d shows 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 3 , respectively, while the density measured by the experiment was 0.81 g/cm 3 [32].This was because C 20 was applied in the MD model to uniformly represent the C 20+ 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 by analogy.As shown in Figure 3b,e, the heavier hydrocarbon components generally have higher adsorption density in the first adsorption layer.There is no obvious peak density of C 4 in the first adsorption layer in the weakly hydrophilic system II because the pore wall modified by methyl groups of the weakly hydrophilic system 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 4 adsorbed near the pore wall, which is 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).
Molecules 2024, 29, x FOR PEER REVIEW 6 of 25 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 3 ns 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 about 5 Å. Figure 3a,d shows that the average bulk density of shale oil calculated by MD simulation from  = −11 Å to  = 11 Å was 0.76 and 0.758 g/cm 3 , respectively, while the density measured by the experiment was 0.81 g/cm 3 [32].This was because C20 was applied in the MD model to uniformly represent the C20+ 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 by analogy.As shown in Figure 3b,e, the heavier hydrocarbon components generally have higher adsorption density in the first adsorption layer.There is no obvious peak density of C4 in the first adsorption layer in the weakly hydrophilic system II because the pore wall modified by methyl groups of the weakly hydrophilic system 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 C4 adsorbed near the pore wall, which is 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, it is less clear whether the reduction in the number of different hydrocarbon components in the first adsorption layer is achieved by slipping to the outside of the pores or by desorption from the pore walls into the bulk phase fluid.To evaluate the transport characteristics of different hydrocarbon components of adsorbed oil in the pores, the migration trajectories of C 4 , C 8 , and C 20 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, it is less clear whether the reduction in the number of different hydrocarbon components in the first adsorption layer is achieved by slipping to the outside of the pores or by desorption from the pore walls into the bulk phase fluid.To evaluate the transport characteristics of different hydrocarbon components of adsorbed oil in the pores, the migration trajectories of C4, C8, and C20 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 () 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, was calculated [38].
where  () is the position of molecule i at time t,  (0) is the initial position of molecule i, and the signal 〈…〈 represents an ensemble average.For the hydrophilic system, the  in the X direction of the hydrocarbons is much larger than that in the Z direction (Figure 5a,c), 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 in SI II decreases, the  in the X direction of the first adsorption layer of shale oil decreases and is similar to the  in the Z direction (Figure 5b,d), which means that the slip effect along the wall surface weakens and the 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 (Figure 3b,e), leading to a slightly enhanced diffusion 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, was calculated [38].
where r i (t) is the position of molecule i at time t, r i (0) is the initial position of molecule i, and the signal ⟨. ..⟩ represents an ensemble average.For the hydrophilic system, the MSD in the X direction of the hydrocarbons is much larger than that in the Z direction (Figure 5a,c), 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 in SI II 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 5b,d), which means that the slip effect along the wall surface weakens and the 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 (Figure 3b,e), leading to a 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 5c,d).Both in the strongly and weakly hydrophilic systems, the heavier hydrocarbon component has a smaller MSD in the X and Z 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 SI I is more significant is because the adsorption effect of C 20 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.
of light components, but there is no significant change in the  in the Z direction of the medium and heavy components (Figure 5c,d).Both in the strongly and weakly hydrophilic systems, the heavier hydrocarbon component has a smaller  in the X and Z 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 SI I is more significant is because the adsorption effect of C20 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.

Microscopic Advancement Mechanism of Imbibition Front
The driving force for the 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 between water molecules and the pore wall (Figure S2a,b).However, the approximately consistent change characteristics of Figure S2a,b indicate that the 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 [39].Furthermore, as shown in Figure 2a,b, the oil phase has a larger  in more hydrophilic pores; however, the slip of shale oil on the wall surface is significantly reduced

Microscopic Advancement Mechanism of Imbibition Front
The driving force for the 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 between water molecules and the pore wall (Figure S2a,b).However, the approximately consistent change characteristics of Figure S2a,b indicate that the 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 [39].Furthermore, as shown in Figure 2a,b, the oil phase has a larger L S 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 L S 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 [40] and visualization results from our MD simulation, 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.
in the weakly hydrophilic system.Therefore, the wettability limitation of the  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 [40] and visualization results from our MD simulation, 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.As shown in Figure 6c, 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 As shown in Figure 6c, 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 20 H 42 ) (M abo ) in SI 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 20 H 42 in Figure 4.
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,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 (The positional coordinates over time of "interfacial water molecule" in the X and Z directions in weakly hydrophilic systems also illustrate this issue, as shown in Figure 6f,g).
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 displacement efficiency (DE) [41]: where N O represents the number of oil molecules originally in the pore; N t represents the number of oil molecules in the pore at time t.
As shown in Figure 7d, we observe that the DE increases approximately linearly in the first stage I (At 0-1621 ps).However, the increase in DE 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 (inset in Figure 7d).The DE in stage I amounts to 69.23%, with a DE per unit time of 0.043% ps −1 .Moreover, the DE in stage II amounts to 30.77%, with a DE 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 DE 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 SI I and SI II at 500 ps in Figure 7a,b show significant fluctuations at the oil-water interface.However, the velocity fluctuations of water molecules at the remaining pore positions, except for the oil-water interface region in the X direction, are 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). ecules".
In addition, the velocity distribution of water molecules in the X direction in both SI I and SI II at 500 ps in Figure 7a,b show significant fluctuations at the oil-water interface.However, the velocity fluctuations of water molecules at the remaining pore positions, except for the oil-water interface region in the X direction, are 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 SI I and SI II were about 0.275 nm away from the hydroxyl groups on the wall, which was close to the distance of the hydrogen bonds [42].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 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 SI I and SI II were about 0.275 nm away from the hydroxyl groups on the wall, which was close to the distance of the hydrogen bonds [42].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, i.e., water molecules spread forward along the pore walls.By comparing the g(r) values in Figure 8a,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.
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, i.e., water molecules spread forward along the pore walls.By comparing the g(r) values in Figure 8a,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.

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 [43].To simulate the forced imbibition process when the pressure difference (∆, unlike ∆ in Equation (3), represents only the difference between shut-in pressure and original matrix pore pressure, the same as below) 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 two forced imbibition systems corresponding to systems I and II are named forced imbibition system I (FI I) and forced imbibition system II (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)

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 [43].To simulate the forced imbibition process when the pressure difference (∆P, unlike ∆P in Equation ( 3), represents only the difference between shut-in pressure and original matrix pore pressure, the same as below) 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 two forced imbibition systems corresponding to systems I and II are named forced imbibition system I (FI I) and forced imbibition system II (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).
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  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 change is relatively large, and the diffusion ability of the medium component in the Z direction changes relatively little (Figure 10a,b,d,e).The  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 irregular changes in adsorption diffusion properties but does not effectively strip the adsorbed layer of shale oil (especially the heavy component).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 change is relatively large, and the diffusion ability of the medium component in the Z direction changes relatively little (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 irregular changes 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 in the pressure difference, the degree of curvature of the meniscus at the imbibition front tends to weaken (i.e., the value of θ d increases with the increase in shut-in pressure), and it is especially obvious at 20 MPa.In the weakly hydrophilic system, 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 a pressure difference [44].
Further observation of the imbibition front of the simulated snapshots of forced imbibition in Figure 9a,b shows that with the increase in the pressure difference, the degree of curvature of the meniscus at the imbibition front tends to weaken (i.e., the value of  increases with the increase in shut-in pressure), and it is especially obvious at 20 MPa.In the weakly hydrophilic system, 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 a pressure difference [44].Because the uneven distribution of mixed crude oil results in an asynchronous advancement of the imbibition front and a 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 the velocity distribution of shale oil using Equation (3), as shown in Figure 11a.Considering the asynchronous advancement of the imbibition front and the accuracy of the results, both the slip velocity ( ) and peak velocity ( ) are obtained by arithmetic averaging, as shown in Equations ( 6) and (7).
where  ,  , and  are the slip velocities of the shale oil on the lower and upper walls, and their arithmetic mean, respectively;  ,  , and  are the Because the uneven distribution of mixed crude oil results in an asynchronous advancement of the imbibition front and a 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 the velocity distribution of shale oil using Equation (3), as shown in Figure 11a.Considering the asynchronous advancement of the imbibition front and the accuracy of the results, both the slip velocity (v slip ) and peak velocity (v peak ) are obtained by arithmetic averaging, as shown in Equations ( 6) and (7).
where v lower slip , v upper slip , and v slip are the slip velocities of the shale oil on the lower and upper walls, and their arithmetic mean, respectively; v lower peak , v upper peak , and v peak are the peak velocities of shale oil on the lower and upper sides of the center of the pore, and their arithmetic mean, respectively.The fitting results for systems I and II are presented in Figures 11b and 11c, respectively.respectively, show that the ratio of the weakly hydrophilic system is greater than that of the strongly hydrophilic system.This finding suggests that the difference between the  and  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.Linear fitting again of v slip and v peak obtained by Equations ( 6) and ( 7) can obtain the slopes K peak and K slip , respectively representing the change rate of peak velocity and slip velocity with pressure difference.Both v slip and v peak increase faster for the strongly hydrophilic system than for the weakly hydrophilic system, 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 an increase in the dynamic contact angle (i.e., wetting hysteresis).Further comparison of the ratios of K peak and K slip (i.e.,

K peak K slip
) for the different imbibition systems, respectively, show that the ratio K peak K slip of the weakly hydrophilic system is greater than that of the strongly hydrophilic system.This finding suggests that the difference between the v peak and v slip 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 [45], when forced imbibition occurs after shut-in, the equilibrium of forces on the fluid in the confined capillary can be expressed as: where F c is the capillary force; F P 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; and 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 [46].Secondly, for pores with extremely small radii, f dominates and the F i can be ignored [47].Then Equation ( 8) can be simplified to the following: Assuming that the width between two walls, length in the Y direction, and length in the X direction of the slit-pore are w, l, and L, respectively, it is possible to convert displacement force F P into wl∆P.According to the Young-Laplace Equation [48], the F c in a single slit-pore is the following: where θ is the contact angle of the oil-water-wall three-phase system, which is assumed to be equal to the test results of 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: 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 the following [49]: where v is the imbibition velocity, and µ w and µ o are the viscosity of water and oil, respectively.After bringing Equations ( 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 the following: After integrating Equation ( 13) under the initial condition of x = 0 and t = 0 (dt are the integration variable), we get the following: Additionally, it should be noted that the shale oil model used in this study does not account for components such as C 20+ and asphaltenes, implying that its viscosity could be lower under reservoir conditions.As previously mentioned, the oil DE 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, i.e., the viscosity of water under the same conditions obtained from the NIST database.Consequently, Equation ( 14) can be simplified as follows: Then the imbibition velocity at pressure difference is The pressure differences are equal to 0, 5, 10, and 20 MPa, respectively, in this study.When the pressure difference is equal to 0, Equation ( 16) is converted to v = σwcos θ/12µL, indicating that the spontaneous imbibition velocity is constant during spontaneous imbibition under reservoir conditions.This may be the reason why the oil DE changes approximately linearly under reservoir conditions before the imbibition front reaches the outlet.Some parameters of Equation ( 16) are: The σ under reservoir conditions obtained by MD simulation is 55.44 mN/m (Figure S3) (The results indicate a good match between the simulation and experimental value of 53.2 mN/m [50]); The static contact angles (θ S ) for the strongly and weakly hydrophilic walls are 53.9 • and 78.97 • , respectively (Figure S5d); The L and w are 15 and 6 nm (see the Section 3 and Figure S4).
The imbibition velocity when the pressure difference is equal to 0 is called the spontaneous imbibition velocity (v SI ), and the imbibition velocity when the pressure difference is equal to 5, 10, and 20 MPa is called the forced imbibition velocity (v FI ).Assuming that, when capillary imbibition occurs, the θ does not change with other factors such as the pressure difference, the ratio (R) of forced imbibition to spontaneous imbibition velocity at the pressure difference is as follows: From Equation ( 17), as assumed earlier in the subsection, the three-phase contact angle θ is not affected by the pressure difference and is always equal to the θ S .Therefore, the R of the weakly hydrophilic system should be greater than that of the strongly hydrophilic system at the same pressure differences.However, the calculated R 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 difference leads to the occurrence of the wetting hysteresis phenomenon, i.e., θ changes under the influence of the pressure difference.The imbibition velocity at different pressure differences obtained by substituting the value of these parameters mentioned above into Equation ( 16) is renamed as the theoretical imbibition velocity (v th ) under different pressure differences.Meanwhile, we rename the imbibition velocity obtained through the MD simulation in Table 1 as the MD simulation velocity (v MD ).The calculation method of the imbibition velocity obtained through MD simulation involves 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 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 in capillary pressure due to wetting hysteresis.Consequently, excessively high imbibition pressure impedes the maximization of economic benefits and reduces the underutilization of natural power.
Molecules 2024, 29, x FOR PEER REVIEW 18 of 25 both  and  gradually increase, while ∆ 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, ∆ 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 in capillary pressure due to wetting hysteresis.Consequently, excessively high imbibition pressure impedes the maximization of economic benefits and reduces the underutilization of natural power.

Model System
All molecular models were constructed through Materials Studio software (version: Materials Studio 2020) developed by Accelrys Company in the San Diego, CA, USA.The molecular model constructed in this article mainly includes solid quartz walls, nanopores composed of quartz walls, and fluid components (shale oil and water).
China's nonmarine shale is typically composed of variable amounts of detrital minerals such as quartz, feldspar, and calcite with content generally higher than 40%, and clay minerals such as smectite, illite, and kaolinite [51].In addition, the brittle mineral content affects the microcrack development, oil-bearing capacity, and fracturing stimulation techniques.Higher content of fragile minerals such as quartz and feldspar, and lower content of clay minerals such as kaolinite and smectite, enhance rock brittleness.High brittleness means that intrinsic fractures and induced fractures are prone to appear in the shale under the action of external forces, which would facilitate shale oil recovery [52,53].Therefore, quartz is one of the key matrix minerals in the shale oil extraction process.
The research object of this study is the shale oil reservoirs of the Lucaogou Formation in the Jimsar Sag, Xinjiang oilfield.The reservoir has the geological characteristics of integrated source and storage, and a dispersed "sweet spot", with an average reservoir temperature of 80 °C [54,55] and an original formation pressure of about 40 MPa [56,57].The Luchaogou Formation (P2l) consists of upper (P2l2) and lower (P2l1) sections, each with its own "sweet spot" enrichment sections.Scholars [58,59] have conducted mineral components analysis on shale core samples from the Jimsar reservoir by using a device of Xray 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  China's nonmarine shale is typically composed of variable amounts of detrital minerals such as quartz, feldspar, and calcite with content generally higher than 40%, and clay minerals such as smectite, illite, and kaolinite [51].In addition, the brittle mineral content affects the microcrack development, oil-bearing capacity, and fracturing stimulation techniques.Higher content of fragile minerals such as quartz and feldspar, and lower content of clay minerals such as kaolinite and smectite, enhance rock brittleness.High brittleness means that intrinsic fractures and induced fractures are prone to appear in the shale under the action of external forces, which would facilitate shale oil recovery [52,53].Therefore, quartz is one of the key matrix minerals in the shale oil extraction process.
The research object of this study is the shale oil reservoirs of the Lucaogou Formation in the Jimsar Sag, Xinjiang oilfield.The reservoir has the geological characteristics of integrated source and storage, and a dispersed "sweet spot", with an average reservoir temperature of 80 • C [54,55] and an original formation pressure of about 40 MPa [56,57].The Luchaogou Formation (P 2 l) consists of upper (P 2 l 2 ) and lower (P 2 l 1 ) sections, each with its own "sweet spot" enrichment sections.Scholars [58,59] 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 and is generally hydrophilic [41,[60][61][62].Furthermore, it is mentioned in physics of petroleum reservoirs that spontaneous imbibition can only occur in hydrophilic capillaries.Therefore, considering that minerals themselves were not one of the research variables, in summary, a single mineral quartz is used in this work to represent shale reservoir matrix pores.The initial quartz unit cell SiO 2 Quartz was taken from the structure database of Materials Studio software and then was cleaved along the (010) 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 [63].To characterize the different hydrophilicity of shale reservoirs, silica was modified concerning existing research methods.The strongly 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 010 crystal surface (surface in contact with oil), with a surface hydroxyl density of 9.6 nm −2 , which was consistent with the results of the crystal chemistry calculations (5.9-18.8nm −2 ) [64].Then, hydroxyl hydrogens of the 010 crystal surface were substituted with a certain amount of hydrophobic group methyl (-CH 3 ) to obtain the weakly hydrophilic wall surface (named W B ).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 −2 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 (i.e., hydrophobic) pores cannot occur effectively, only hydrophilic matrix pores of the shale 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 13.The nanopore representing a slit-shaped quartz nanopore in the shale matrix consisted of two quartz walls with the same modification, as shown in Figure S4.Zones A and B represent microfracture or macropore to which the pore is connected at the water-injection end and microfracture or macropore to which the pore is connected at the non-water-injection end, respectively.Furthermore, previous studies have shown that the Jimsar shale reservoir develops abundant nanopores [65][66][67].The Young-Laplace equation [48] shows capillaries with smaller radii have more significant imbibition effects.Secondly, MD simulations are subject to limited computing resources.Therefore, comprehensively considering the computational resources and the significant imbibition effect, the vertical distance (w) between their wall surfaces in the Z direction is set to 60 Å, and the two walls of the pore were placed parallel to the X-Y plane.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 (n-C 4 H 10 ), n-Octane (n-C 8 H 18 ), and n-Eicosane (n-C 20 H 42 ) 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 4 H 10 ), n-Octane (n-C 8 H 18 ), and n-Eicosane (n-C 20 H 42 ) respectively represent the light, medium, and heavy components in Jimsar 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 injected water in the microfracture or macropore to which the pore is connected, 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 3 ns equilibrium molecular dynamics (EMD) at reservoir conditions (353.15K and 40 MPa), thereby obtaining a reasonable shale oil density.During the imbibition simulation, He plates, which are parallel to the Z-Y plane, 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.Furthermore, plates 1 and 2 can only be moved parallel to direction X.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 40 MPa, while the pressure exerted on He plate 1 representing the shut-in pressure (P s ) is greater than or equal to P m .When the pressure on He plate 1 is equal to P m , water is spontaneously imbibed into the pore.When the pressure applied to He plate 1 exceeds P m , it is the process of forced imbibition oil displacement.The simulation time of imbibition oil displacement is 1 ns.process was controlled by controlling the pressure attached to the He plate.Furthermore, plates 1 and 2 can only be moved parallel to direction X.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 ( ) of 40 MPa, while the pressure exerted on He plate 1 representing the shut-in pressure ( ) is greater than or equal to  .When the pressure on He plate 1 is equal to  , water is spontaneously imbibed into the pore.When the pressure applied to He plate 1 exceeds  , it is the process of forced imbibition oil displacement.The simulation time of imbibition oil displacement is 1 ns.

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. [68] for hydrates and multi-component minerals and their interaction with fluid interfaces and has been successfully applied to hydroxylated silica surface-fluid systems [69][70][71][72].Thus, 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 [73].The OPLS-AA force field [74] is used to describe n-hydrocarbons in shale oil.The OPLS force field can accurately describe hydrocarbon molecules' thermodynamic properties and structural characteristics [75].The SPC/E potential energy model [76]

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. [68] for hydrates and multi-component minerals and their interaction with fluid interfaces and has been successfully applied to hydroxylated silica surface-fluid systems [69][70][71][72].Thus, 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 [73].The OPLS-AA force field [74] is used to describe n-hydrocarbons in shale oil.The OPLS force field can accurately describe hydrocarbon molecules' thermodynamic properties and structural characteristics [75].The SPC/E potential energy model [76] 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 ClayFF force field and OPLS-AA force field [77][78][79][80].The He plate and quartz wall 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 (version: 28 Mar 2023) to perform MD simulations [81], with a time step of 1fs.The simulation results were visualized using Visual Molecular Dynamics (VMD) software (version: 1.9.4) [82].Intermolecular non-bonded interactions (E nonbonded ) can be described by the standard 12-6 Lennard-Jones potential (F LJ ) and Coulomb potential, respectively (E coulomb ) [70,83].
 + 1 4πε 0 ε r q i q j 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; ε 0 and ε r represent the dielectric constant of a vacuum and relative permittivity of the medium, respectively.And q i and q j represent the charge of atoms i and j, respectively.The cutoff radius was set to be 10 Å.All MD simulations were performed in the NVT ensemble.The interactions between different atoms are calculated using the Lorentz-Berthelot mixing criterion [84]:

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 on imbibition 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 improving imbibition oil displacement methods during the shut-in period, optimizing shut-in time and operation after the shale reservoir has been hydraulically fractured, and eventually enhancing production in shale oil reservoirs.

Supplementary Materials:
The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/molecules29051112/s1, Figure S1: The interaction energy between oil and wall during structural optimization of the oil-pore systems of imbibition systems I and II; Figure S2 S1: The average mass fraction of light, heavy, and medium hydrocarbon components of the Jimsar shale oil; Table S2: The mole fractions and number of each component of the shale oil model; Table S3: Atomic charge and detailed interaction parameter of SiO 2 , oil, and water.References [59,63,85,86] are cited in the Supplementary Materials.

Figure 1 .
Figure 1.(a) Simulation snapshots of SI I and SI II, and (b) corresponding dynamic contact angles at several different moments.In the snapshots of the imbibition systems in (a), the blue zone represents water molecules, the green zone represents the SiO2 wall surface, and the gray zone represents hydrocarbon molecules, respectively, the same as below.The inset in (b) shows a schematic diagram of the  on the upper pore wall and  on the lower pore wall.

Figure 1 .
Figure 1.(a) Simulation snapshots of SI I and SI II, and (b) corresponding dynamic contact angles at several different moments.In the snapshots of the imbibition systems in (a), the blue zone represents water molecules, the green zone represents the SiO 2 wall surface, and the gray zone represents hydrocarbon molecules, respectively, the same as below.The inset in (b) shows a schematic diagram of the θ u on the upper pore wall and θ l on the lower pore wall.
dv dz z surf is the velocity gradient of the fluid at the fluid-solid interface.

Figure 2 .
Figure 2. Velocity distribution of shale oil in the Z direction and L S obtained by fitting velocity distribution in the (a) SI I and (b) SI II, respectively.(c) 2-D density distribution of shale oil at different moments during spontaneous imbibition in SI I and SI II.The purple zone in (c) represents the background of the 2-D density distribution of shale oil, the long white bar represents the quartz wall surface, and the rest of the zone is the 2-D density distribution of shale oil, the same as below.

Figure 3 .
Figure 3.The mass density distribution of shale oil in the pores of (a) SI I and (d) SI II after 3 ns of structural optimization, respectively.The number density distribution of different hydrocarbon components in the pores of (b) SI I and (e) SI II after 3 ns of structural optimization, respectively.The number of changes of different hydrocarbon components in the first adsorption layer with time during spontaneous imbibition in (c) SI I and (f) SI II, respectively.In (a,b,d,e), the long purplish red bar represents the first adsorption layer.

Figure 3 .
Figure 3.The mass density distribution of shale oil in the pores of (a) SI I and (d) SI II after 3 ns of structural optimization, respectively.The number density distribution of different hydrocarbon components in the pores of (b) SI I and (e) SI II after 3 ns of structural optimization, respectively.The number of changes of different hydrocarbon components in the first adsorption layer with time during spontaneous imbibition in (c) SI I and (f) SI II, respectively.In (a,b,d,e), the long purplish red bar represents the first adsorption layer.

Figure 4 .
Figure 4. Migration trajectories of C4, C8, and C20 in the first adsorption layer after structural optimization during spontaneous imbibition.The red molecules in the blue circle are labeled as different hydrocarbon molecules in the first adsorption layer after structural optimization, and the yellow, tan, and purple molecules represent C4, C8, and C20, respectively, the same as below.

Figure 4 .
Figure 4. Migration trajectories of C 4 , C 8 , and C 20 in the first adsorption layer after structural optimization during spontaneous imbibition.The red molecules in the blue circle are labeled as different hydrocarbon molecules in the first adsorption layer after structural optimization, and the yellow, tan, and purple molecules represent C 4 , C 8 , and C 20 , respectively, the same as below.

Figure 5 .
Figure 5.  in (a) X and (b) Z directions of different hydrocarbon components in the first adsorption layer after structural optimization in SI I during spontaneous imbibition, respectively. in (c) X and (d) Z directions of different hydrocarbon components after structural optimization in the first adsorption layer in SI II during spontaneous imbibition, respectively.

Figure 5 .
Figure 5. MSD in (a) X and (b) Z directions of different hydrocarbon components in the first adsorption layer after structural optimization in SI I during spontaneous imbibition, respectively.MSD in (c) X and (d) Z directions of different hydrocarbon components after structural optimization in the first adsorption layer in SI II during spontaneous imbibition, respectively.

Figure 6 .
Figure 6.(a) Migration trajectories of "interfacial water molecules" and "adsorbed oil molecules" in SI I. (b) Migration trajectories of "interfacial water molecules" in SI II; (c) Schematic of the advancing imbibition front.Position coordinates of two "interfacial water molecules" (Miw1 and Miw2) and one "adsorbed oil molecule" (Mabo) in the (d) X and (e) Z directions in SI I at several different moments, respectively.Position coordinates of one "interfacial water molecule" in the (f) X and (g) Z directions in SI II at several different moments, respectively.

Figure 6 .
Figure 6.(a) Migration trajectories of "interfacial water molecules" and "adsorbed oil molecules" in SI I. (b) Migration trajectories of "interfacial water molecules" in SI II; (c) Schematic of the advancing imbibition front.Position coordinates of two "interfacial water molecules" (M iw1 and M iw2 ) and one "adsorbed oil molecule" (M abo ) in the (d) X and (e) Z directions in SI I at several different moments, respectively.Position coordinates of one "interfacial water molecule" in the (f) X and (g) Z directions in SI II at several different moments, respectively.

Figure 7 .
Figure 7. Velocity distribution of water molecules in the X direction in (a) SI I and (b) SI II at 500 ps, respectively.(c) Schematic diagram of the migration trajectories of water molecules outside the pore.(d) Imbibition  of SI I.The purplish red bar represents the oil-water interface region in (a,b).Red, purple, dark blue, and light blue balls in (c) represent water molecules outside the pore successively adsorbed to the pore wall, respectively, and the black balls represent hydrocarbon molecules.The inset in the red circle in (d) means the imbibition front reaches the right end of the nanopore on the upper walls.

Figure 7 .
Figure 7. Velocity distribution of water molecules in the X direction in (a) SI I and (b) SI II at 500 ps, respectively.(c) Schematic diagram of the migration trajectories of water molecules outside the pore.(d) Imbibition DE of SI I.The purplish red bar represents the oil-water interface region in (a,b).Red, purple, dark blue, and light blue balls in (c) represent water molecules outside the pore successively adsorbed to the pore wall, respectively, and the black balls represent hydrocarbon molecules.The inset in the red circle in (d) means the imbibition front reaches the right end of the nanopore on the upper walls.

Figure 8 .
Figure 8. Radial distribution function g(r) of O (hydroxyl) and O (H2O) atoms in the (a) SI I and (b) SI II, respectively.g(r) of O (H2O) and C (hydrocarbons) atoms in the (c) SI I and (d) SI II.

Figure 8 .
Figure 8. Radial distribution function g(r) of O (hydroxyl) and O (H 2 O) atoms in the (a) SI I and (b) SI II, respectively.g(r) of O (H 2 O) and C (hydrocarbons) atoms in the (c) SI I and (d) SI II.

Figure 9 .
Figure 9. Snapshots of (a) FI I and (b) FI II under different pressure differences at several moments.(c) 2-D density distribution of shale oil at 1000 ps in FI I and FI II.The purple zone in (c) represents the background of the 2-D density distribution of shale oil, the long white bar represents the quartz wall surface, and the rest of the zone is the 2-D density distribution of shale oil.

Figure 9 .
Figure 9. Snapshots of (a) FI I and (b) FI II under different pressure differences at several moments.(c) 2-D density distribution of shale oil at 1000 ps in FI I and FI II.The purple zone in (c) represents the background of the 2-D density distribution of shale oil, the long white bar represents the quartz wall surface, and the rest of the zone is the 2-D density distribution of shale oil.

Figure 10 .
Figure 10. in the Z direction for (a) C4, (b) C8, and (c) C20 in the first adsorption layer after structural optimization in FI I during forced imbibition, respectively. in the Z direction for (d) C4, (e) C8, and (f) C20 in the first adsorption layer after structural optimization in FI II during forced imbibition, respectively.

Figure 10 .
Figure 10.MSD in the Z direction for (a) C 4 , (b) C 8 , and (c) C 20 in the first adsorption layer after structural optimization in FI I during forced imbibition, respectively.MSD in the Z direction for (d) C 4 , (e) C 8 , and (f) C 20 in the first adsorption layer after structural optimization in FI II during forced imbibition, respectively.

Figure 11 .
Figure 11.(a) Velocity distribution of shale oil and schematic of slip velocity ( and  ) and peak velocity ( and  ) obtained by fitting velocity distribution of shale oil.(b)  and  of shale oil at different pressure differences, and slope obtained by linearly fitting  and  of shale oil in FI I. (c)  and  of shale oil at different pressure differences, and slope obtained by linearly fitting  and  of shale oil in FI II.

Figure 11 .
Figure 11.(a) Velocity distribution of shale oil and schematic of slip velocity (v lower slip and v upper slip ) and peak velocity (v lower peak and v upper peak ) obtained by fitting velocity distribution of shale oil.(b) v slip and v peak of shale oil at different pressure differences, and slope obtained by linearly fitting v slip and v peak of shale oil in FI I. (c) v slip and v peak of shale oil at different pressure differences, and slope obtained by linearly fitting v slip and v peak of shale oil in FI II.

Figure 12 .
Figure 12. v MD , v th , and ∆v in (a) strongly hydrophilic and (b) weakly hydrophilic systems, respectively.

3 .
Models and Methodology 3.1.Model System All molecular models were constructed through Materials Studio software (version: Materials Studio 2020) developed by Accelrys Company in the San Diego, CA, USA.The molecular model constructed in this article mainly includes solid quartz walls, nanopores composed of quartz walls, and fluid components (shale oil and water).

Figure 13 .
Figure 13.Snapshot of the initial configuration of the imbibition model.
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 ClayFF force field and OPLS-AA force field [77-80].The He plate and quartz wall 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 (version: 28 Mar 2023) to perform MD simulations[81], with a time step of 1fs.The simulation results were visualized using Visual Molecular Dynamics (VMD) software (version: 1.9.4)[82].Intermolecular non-bonded interactions ( ) can be described by the standard 12-6 Lennard-Jones potential ( ) and Coulomb potential, respectively ( )[70,83].

Figure 13 .
Figure 13.Snapshot of the initial configuration of the imbibition model.
: (a) The interaction energy between water molecules and wall, and (b) Hy-drogen bond number as a function of time t in spontaneous imbibition systems I-II, respectively; Figure S3: Molecular model for calculating the interfacial tension of pure oil and water systems; Figure S4: Quartz wall sizes and different modified wall surfaces with different wettability; Figure S5: (a) Initial and (b) equilibrium configuration for static contact angle testing, respectively.(c) Schematic diagram of static contact angle calculation method.(d) Static wetting angle of W A and W B ; Figure S6: The mass fraction of hydrocarbons with different numbers of carbon atoms of (a) JHW05815 Oil Sample, (b) JHW07121 Oil Sample, and (c) J41 Oil Sample; Table

Table 1 .
v MD and R under different pressure differences in imbibition systems I and II.