Pores in tight oil reservoirs are mainly nano-scale, so the fluid occurrence and flow mechanism are different from conventional oil reservoirs ^[1]. At present, fracturing stimulation and water injection huff-puff are the primary methods for developing tight oil reservoirs. In the development process, influenced by capillary force, water (fracturing fluid) in fractures flows into matrix and displaces oil ^[2,3,4,5].

The imbibition parameters of nano-scale pores are mainly affected by nanoconfinement effect. Wu et al. ^[6] studied the nanoconfinement effect in pores on the viscosity and slippage effect of water and obtained a Poiseuille equation under the conditions of nano-scale single-pore and single-phase. Supple et al. ^[7] carried out an imbibition experiment on liquid flooding gas in nano- scale pores and believed that the imbibition does not follow the LWR (Lucas-Washburn-Rideal) equation ^[8], and the imbibition distance is linear with time. Feng et al. ^[9] discussed the change of gas-water interfacial tension in nano-scale pores, but did not apply it to a nano-scale imbibition process. Fernández-Toledano et al. ^[10] validated the applicability of the molecular kinetic theory for nano- scale pores ^[11] in a solid-gas-liquid imbibition system. Kelly et al. ^[12] studied the imbibition problem of gas-liquid in nano-scale pores and found that due to the influence of liquid film and fluid viscosity changes, the imbibition distance cannot be calculated by the traditional LWR equation. Gruener et al. ^[13] studied the influence of nanoconfinement on gas-liquid imbibition in cores and realized that the continuity of the meniscus is affected by the length of the pore. Liu et al. ^[14] found a linear relationship between recovery and time in initial imbibition when studying oil-water imbibition in tight cores. In conclusion, less studies have carried out on oil-water imbibition in nano-scale pores, and the mechanism of dynamic contact angle on oil-water imbibition in nano-scale pores is not clear.

In the present study, considering the effects of dynamic contact angle, nanoconfinement (including multilayers sticking and slippage), inertia, and inlet end, an equation for oil-water imbibition in nano-scale pores is established, the relationship between the friction coefficient of solid-oil-water three-phase contact line and fluid viscosity on the contact is derived. In addition, based on the capillary bundle model and the lognormal distribution theory, the imbibition model for tight cores is upgraded. Finally, the adaptability of the model is validated by imbibition simulation to pores and tight cores, and key parameters affecting the imbibition process are analyzed.

Based on the LWR theory, the oil-water imbibition equation in conventional pores (marked a LWR model) ^[8] is as follows:

Considering nanoconfinement effect, inlet end effect, and inertia effect, the imbibition equation in nanopores (marked a ND model) is expressed as ^[15]:

The ND model ignores the effect of dynamic contact angle (that is, it is assumed that contact angle does not change with imbibition velocity). However, with friction effect on contact line, contact angle varies with imbibition velocity, and the process is dynamic. If considering the effect of dynamic contact angle, the imbibition equation in nanopores (marked a CD model) is expressed as:

The inertial force in the equation is:

The capillary force calculated by the equilibrium contact angle is ^[15]:

Based on the Poiseuille equation considering slippage, the viscous force can be expressed as ^[6]:

The slip length can be expressed as ^{[6, 16]}:

The inlet end effect can be expressed as ^[17]:

The gravity is ^[15]:

The friction force along the three-phase contact line is ^[15]:

In Eq. (12), the friction coefficient of the three-phase contact line characterizes the strength of the dynamic contact angle effect.

According to molecular kinetic theory, the relationship between dynamic contact angle and velocity is ^[18,19]:

Considering that the molecules of solid phase, wetting phase and non-wetting phase exist in a mixed form in the three-phase contact zone, a weighting method is proposed to calculate the wetting activation free energy:

In the equation, the weight coefficient is the ratio of the viscosity of single-phase fluid to total viscosity.

The activation free energy of viscosity effect is characterized by the Eyring theory ^[20]:

The activation free energy of rock-fluid interaction is expressed by adhesion work ^[21,22]:

From Eqs. (14)-(19), it can be seen that the friction coefficient of the three-phase contact line is:

Considering multilayer sticking effect, the viscosity of the solid-liquid interface area is:

Combining Eqs. (18) and (19), the fluid viscosity in the interface area is:

Combining Eqs. (20), (22), and (23), the relationship between the friction coefficient of the contact line and the fluid viscosity in the interface area can be obtained:

In addition, the effective fluid viscosity can be calculated by the area weight method ^[6].

Wang et al. ^[23] studied the effects of molecular interaction, dynamic contact angle and inlet end on water imbibition dynamics in nanoporous media and obtained relatively reliable imbibition data (simulation data validated by Lattice Boltzmann). To discuss the degree of influence of multilayer sticking, slippage, inlet end, and inertia on imbibition dynamics, the imbibition parameters in this reference are used as the basic parameters of the CD model in this paper (Table 1), and the imbibition dynamics data is used as validation data to simulate and analyze the effects of multilayers sticking, slippage, inlet end, and inertia on the result of imbibition.

The effective interfacial tension is calculated as follows ^[24]:

The thickness of the interface area is the thickness of two layers of fluid molecules ^[6], and the density of the wetting phase in the interface area is calculated as follows ^[25]:

Assuming that the ratio of the density of the wetting phase in the interface area to the conventional density of the fluid and the ratio of the density of the non-wetting phase in the interface area to the conventional density of the fluid are the same, the density of the non-wetting phase in the interface area is also obtained.

The effective density is calculated using the weight method:

Fig. 1 shows the comparison curves between simulation results and validation data when the CD model ignores the effects of multilayers sticking, slippage, inlet end, inertia, etc. It can be seen that when inlet end and inertial effects are ignored, the simulation curve coincides with the validation curve, indicating that inlet end and inertia effects have a small impact on imbibition and can be ignored. When multilayers sticking effect is ignored, the simulation curve deviates far from the validation curve, and the error is significant, indicating multilayers sticking effect can’t be ignored.

Calculations show the slip length of water in hydrophilic reservoirs is small, ranging from 0.10 to 0.41 nm; the slip length of oil is related to the properties of the pore wall surface and the composition of the oil. The main component of the oil in this paper is n-octane, and the slip length of the oil is small (0.19 nm). In this case, the impact of the slippage effect on the imbibition dynamics can be ignored. However, when composition changes, the slip length of oil will change. For example, in sandstone with known surface energy of 61.2 mN/m, the slip lengths of n-dodecane and n-tetradecane on the surface are 5.7 nm and 28.0 nm, respectively ^[16]. In this case, the influence of slippage effect on imbibition dynamics is more significant, and it cannot be ignored.

Based on the above analysis, if ignoring the influence of inlet end and inertia effects, we get Eq. (28) from Eqs. (3), (5) - (7), and (12):

To facilitate characterization and calculation, the above equation can be expressed as:

Then the imbibition distance and volume can be obtained as follows:

Considering the influence of connate water saturation and residual oil saturation in the core, the pore volume where fluid can flow is:

Based on the capillary bundles in the core, the pore volume where fluid can flow is:

Combining Eqs. (32) and (33), the total number of capillaries in the core is:

Then the imbibition volume and recovery of the core can be expressed as:

The pore size of the core can be characterized by the lognormal distribution model ^[26]:

The basic data in Table 1 is used to validate the reliability of dynamic contact angle model and pore imbibition model. The specific calculation process is shown in Fig. 2.

Stroberg et al. ^[27] first simulated molecular imbibition dynamics at nanoscale, whose result is used in this study for model validation. First, calculate the fluid viscosity in the interface area, then calculate the friction coefficient along the contact line (where the jump length is used as the fitting parameter), and finally build the relationship between the dynamic contact angle and the velocity by Eq. (13). It can be seen from Fig. 3 that the simulation results of the dynamic contact angle model proposed in this paper and the reference model (molecular dynamic model from Stroberg et al., 27) are in good agreement (the simulated pore radius is 2.5 nm and 4.0 nm, respectively), indicating that the model in this paper can better characterize the dynamic change in contact angle during imbibition process.

Using the basic data in Table 1, the LWR, ND, and CD models are used to simulate the relationship between imbibition time and imbibition distance, and the simulation results are compared with the validation curves (Fig. 4). It can be seen that the simulation curve from the traditional LWR model that does not consider the effects of nanoconfinement, inlet end, inertia and dynamic contact angle is the farthest from the validation curve, and the error is the largest, so it cannot characterize the imbibition process in nano-scale pores; the ND model considers the effects of nanoconfinement, inlet end and inertia, but the effect of dynamic contact angle is not considered. Although the simulation curve is closer to the validation curve than the LWR model, the error is still enormous; based on the ND model, the CD model considers the effect of dynamic contact angle, so the simulation result is better, and the goodness of fit is 0.99. It is clear that dynamic contact angle has the most significant impact on imbibition, followed by nanoconfinement effect (multilayers sticking effect and slippage effect) and inertia effect, the last is inlet end effect which is the most negligible.

Ren et al. ^[28] conducted a study on the influence of permeability and viscosity on the imbibition process in tight cores. The results showed that the recovery increases with the increase of permeability and decreases with the increase of oil phase viscosity. Wang et al. ^[29] studied the effects of temperature and pressure on the imbibition process in tight cores and confirmed that the higher the temperature, the higher the recovery, and there is an optimal pressure range. Kathel et al. ^[30] discussed the influence of wettability on the imbibition process in tight cores and found that surfactants can revert the wettability of the cores, and the spontaneous imbibition recovery is about 68%. A total of 5 experimental cores are selected from the above references, which are core A-5 ^[28], core R4-2 ^[29], core R4-4 ^[29], core 2^[30], and core 7^[30]. They are all tight sandstone cores with basic parameters in Table 2^[28,29,30]. Five sets of data on the imbibition recovery degree with the imbibition time have been obtained through experiments.

According to the data in Table 2, the CD model is used to simulate the relationship between the imbibition recovery degree and the imbibition time, and compared with 5 sets of experimental data (Fig. 5). It can be seen that the simulation results of the CD model are in good agreement with the experimental results from references and can better simulate the oil-water imbibition process in tight cores. However, the CD model ignores the influence of gravity, so when the capillary force is weak, gravity has a certain effect on the imbibition dynamics, resulting in a certain deviation between the simulation curve and the experimental data.

According to the data in Table 1, the ND and CD models are used to simulate how pore radius and interfacial tension affect imbibition dynamics.

Inertia force, contact line friction, the viscous force of the wetting phase (water), and the viscous force of the non-wetting phase (oil) are primary resistances in the imbibition process in tight cores. The ratio of a single resistance to the total of the 4 resistances is defined as a single resistance ratio. Fig. 6 shows the change of single resistance ratios with time in pores with different radii. The following conclusions can be obtained. (1) In initial imbibition, inertia force is the dominant force. Taking 50.0 nm pore as an example, when the imbibition time is less than 0.100 ns, as time increases, the inertial resistance ratio gradually decreases, while the friction resistance ratio gradually increases. When the imbibition time is 0.1-100 ns, the friction resistance ratio is the highest, and the viscous resistance ratio of the oil phase is the second. The friction force along the contact line is dominant, and the friction resistance and the viscous resistance of the oil phase remain unchanged. When the imbibition time is greater than 100 ns, the friction resistance ratio of the contact line and the viscous resistance ratio of the oil phase gradually decrease, while the viscous resistance ratio of the water phase gradually increases. (2) The friction force of the contact line not only acts on the late stage of imbibition, showing a downward trend, but also acts on the initial stage of imbibition, showing an upward trend. (3) As the pore radius decreases, the inertial resistance ratio curve moves to left, and the acting time of inertial force shortens. In the initial stage of imbibition, the friction resistance ratio curve moves to left (corresponding to the simulated curve of the pore radius from 50.0 nm to 20.0 nm and then to 10.0 nm; and limited by the minimum simulation time 0.001 ns in this paper, the simulation curves with pore radii of 5.0 nm and 2.5 nm cannot fully show this process), indicating the effect of the dynamic contact angle increases. In the middle and late stages, the friction resistance ratio curve moves downward, indicating the effect of the dynamic contact angle weakens; the viscous resistance ratio curve of the oil phase moves upward, indicatig the influence of the viscous force of the oil phase increases; the viscous resistance ratio curve of the water phase shifts to left, indicating the initial action time of the viscous force of the water phase advances.

Fig. 7 shows the relationship between imbibition time and imbibition distance and the enhancement factor of imbibition distance (the ratio of the imbibition distance of the CD model to the imbibition distance of the ND model) for different pore radii and on different imbibition models. It can be seen that as the pore radius decreases, the positions of the simulation curve of the CD and ND models drop, and the imbibition distance decreases within the same time. With the extension of the imbibition time, the enhancement factor of the imbibition distance takes a “U” shape, falls in the early stage, stabilizes in the mid-term, and rises in the later period. Because the imbibition distance is negatively correlated with the dynamic contact angle effect, the dynamic contact angle effect rises in the early stage, tends to be stable in the mid-term, and then declines in the later period, which is consistent with the change in the friction resistance ratio of the contact line in Fig. 6b. In addition, as the pore radius decreases, the curve in the initial stage of imbibition moves to left, indicating the influence of the dynamic contact angle increases. The curve in the middle and late stages of imbibition moves upward, showing the impact of the dynamic contact angle is weakened, which is still consistent with the change in the friction resistance ratio of the contact line in Fig. 6b.

Fig. 8 shows the relationship between the imbibition time and the dynamic contact angle for different pore radii. At the beginning of imbibition, the imbibition velocity is approximately zero when the water phase just touches the inner wall of the pore. According to Eq. (13), the contact angle at this time is the equilibrium contact angle (note that the initial time is close to zero, which is much less than the minimum simulation time of 0.001 ns in this paper, so this process is not shown in the Figure). As the imbibition time increases, the imbibition velocity increases, and the contact angle gradually increases to the maximum, and remains stable for a certain period, as shown in Fig. 8. As the imbibition time further increases, the imbibition velocity decreases and the dynamic contact angle gradually decreases, as shown in Fig. 8. According to Eq. (13), it can be predicted that when the imbibition time increases until the imbibition velocity decreases and approaches zero, the contact angle will eventually approach the equilibrium contact angle (and the approaching time is much longer than the maximum simulation time of 1000 ns in this paper).

The greater the oil-water interfacial tension, the greater the imbibition power is, but the adhesion work between fluid and rock will increase too. Fig. 9 shows the relationship between imbibition time and the friction resistance ratio of the contact line under different interfacial tension conditions. It can be seen that as the interface tension increases, the curve moves upward, and the friction resistance ratio of the contact line increases, indicating the dynamic contact angle effect increases.

Fig. 10 is the relationship curve between imbibition time and imbibition distance and the enhancement factor of the imbibition distance under different interfacial tension conditions. It can be seen that as the interfacial tension increases from 0.01 N/m to 0.05 N/m, the positions of the simulation curves of the CD and ND models rise, and the imbibition distances increase within the same imbibition time; but when the interfacial tension increases from 0.05 N/m to 0.10 N/m, the position of the simulation curve of the CD model drops, and the simulation curve of the ND model first rises and then drops. It can be concluded that there is a critical interfacial tension (0.05 to 0.10 N/m in this case). Therefore, if the interfacial tension is too low, it is impossible to obtain a better imbibition effect in tight reservoirs. The best interfacial tension should be from optimization. The results in the late stage of imbibition in Fig. 10c also show that there is a critical interfacial tension.

Fig. 11 shows the relationship curve between imbibition time and dynamic contact angle under different interfacial tension conditions simulated by the CD model. With the extension of the imbibition time, the dynamic contact angle gradually increases from the initial equilibrium contact angle to the maximum and remains stable for a particular period. Then the dynamic contact angle decreases progressively and approaches the equilibrium contact angle. At the same time, the maximum dynamic contact angle increases with the increase in the interfacial tension, that is, as the interfacial tension increases, the change of the contact angle increases, indicating the dynamic contact angle effect increases.

In the imbibition process in nanopores, the effect of dynamic contact angle is the most significant, followed by nanoconfinement (multilayers sticking and slippage), and inertia and inlet end effects are the least.

In early imbibition, the role of inertial force decreases, the role of contact line friction increases, and dynamic contact angle gradually increases from the initial equilibrium to the maximum and remains stable. In late imbibition, the role of contact line friction decreases, the contact angle gradually decreases from the maximum dynamic and approaches the initial equilibrium.

As pore radius decreases, the effect of dynamic contact angle is enhanced in initial imbibition and weakened in middle and late imbibition. As oil-water interfacial tension increases, imbibition power and the effect of dynamic contact angle increase.

There is a critical interfacial tension. Too low interfacial tension cannot bring better imbibition in tight reservoirs. The best interfacial tension should be from optimization.

The model proposed in this paper considers the influences of dynamic contact angle, nanoconfinement, inertia and inlet end, so it can better describe the imbibition process in nanopores and tight reservoirs.

A₁—intermediate variable, Pa•s/m;

A_i—cross-sectional area of ​​the interface area, m²;

A_t—cross-sectional area of ​​the whole flow area, m²;

B—intermediate variable, Pa•s;

C—intermediate variable, N/m;

C_f1—constant related to logarithmic standard deviation of pore radius, m;

C_sw—constant related to the slippage effect of water, m;

d—pore diameter, m;

D_c—core diameter, m;

f₁, f₂—logarithmic standard deviation and logarithmic mean of pore radius, respectively, dimensionless;

f(R)—distribution frequency of pore radius, m^-1;

f(θ_eq)—function related to fluid density and equilibrium contact angle in the solid-liquid interface area, kg/m³;

g—acceleration of gravity, m/s²;

h—Planck's constant, 6.626×10^-34 J•s;

k_B—Boltzmann constant, 1.38×10^-23 J/K;

K—jump frequency, Hz;

l_s—slip length, m;

L—pore length, m;

L_c—core length, m;

n—number of adsorption sites per unit area, m^-2;

N_A—Avogadro constant, 6.02×10²³;

N_m—number of capillaries where fluid can flow;

p—pressure, Pa;

Q—imbibition volume per capillary, m³;

Q_c—core imbibition volume, m³;

R—pore radius, m;

R_ave—average pore radius, m;

R_o—recovery, %;

S—equivalent wall area that allows ​​adjacent molecules enter a pore, m²;

S_wc—connate water saturation, %;

S_or—residual oil saturation, %;

t—imbibition time, s;

T—temperature, K;

v—velocity, m/s;

V—effective volume per fluid molecule, m³;

V_m—pore volume where fluid can flow, m³;

x—imbibition distance, m;

x(R,t)—imbibition distance function, m;

ΔG—activation free energy, J;

α—equivalent wall area ratio composed of solid phase, dimensionless;

γ—interfacial tension, N/m;

δ—average distance between adjacent layers of oil molecules, m;

δ_l—Tolman length, nm;

ζ—friction coefficient of contact line, Pa•s;

θ_d—dynamic contact angle, (°);

θ_eq—equilibrium contact angle, (°);

λ—jump length, m;

μ—viscosity, Pa•s;

ρ—density, kg/m³;

σ₁—surface tension of the oil phase in vacuum, N·m^-1;

τ—tortuosity, dimensionless;

τ_t—relaxation time, s;

φ—pore inclination, (°);

ϕ—porosity, %;

ω—weighting coefficient, dimensionless.

b—fluid parameter value under normal conditions;

ce—calculation based on equilibrium contact angle;

e—effective value of a fluid parameter;

en—inlet end;

fl—contact line friction;

g—gravity;

i—interfacial area;

in—inertia;

l—liquid phase;

max—maximum;

min—minimum;

nw—non-wetting phase;

s—solid phase;

v—viscosity;

w—wetting phase;

wet—wetting.