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.
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 [
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.
where
is the velocity of the fluid at the position of the fluid-solid interface z
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:
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
3 respectively, while the density measured by the experiment was 0.81 g/cm
3 [
33]. 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 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
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 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
4 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
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, 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].
where
is the position of molecule i at time t and
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
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.
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
20H
42) (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
20H
42 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]:
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).
where
,
, and
are the slip velocities of the shale oil on the lower and upper walls, and their arithmetic mean respectively;
,
, and
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 Kpeak and Kslip 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, Kpeak is greater than Kslip, leading to a decrease in the dynamic contact angle (i.e., wetting hysteresis). Further comparison of the ratios of Kpeak and Kslip (i.e., ) for the different imbibition systems, respectively, the ratio 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:
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;
is gravity; Because the matrix pore width of tight shale reservoirs is extremely small and much smaller than the imbibition height, the influence of
can be ignored [
45]. Secondly, for pores with extremely small radii, viscous force
dominates and the inertial force
can be ignored [
46]. Then equation (8) can be simplified to:
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
into wl∆P. According to the Young-Laplace equation [
47], the F
c in a single pore is:
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:
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]:
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:
After integrating equation (13) under the initial condition of x = 0 and t = 0, we get:
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 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:
Then the imbibition velocity at pressure difference is
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:
From equation (17), assuming that the three-phase contact angle
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
under different pressure differences. Meanwhile, we rename the imbibition velocity obtained through MD simulation in
Table 1 as
. 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
and
, as well as the difference (Δv) (∆v=
-
) between the them at different pressure differences. The results show basically that
and
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
and
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.