Phase-Field Simulation of Imbibition for the Matrix-Fracture of Tight Oil Reservoirs Considering Temperature Change

: Injection water temperature is often different from that of the reservoir during water injection development in the tight reservoir. Temperature change causes different ﬂuid properties and oil-water interface properties, which further affects the imbibition process. In this paper, a matrix-fracture non-isothermal oil-water imbibition ﬂow model in tight reservoirs is established and solved by the ﬁnite element method based on the phase-ﬁeld method. The ideal inhomogeneous rock structure model was used to study the inﬂuence of a single factor on the imbibition. The actual rock structure model was used to study the inﬂuence of temperature. The mechanism of temperature inﬂuence in the process of imbibition is studied from the micro-level. It is found that the imbibition of matrix-fracture is a process in which the water enters the matrix along with the small pores, and the oil is driven into the macropores and then into the fractures. Temperature affects the imbibition process by changing the oil-water contact angle, oil-water interfacial tension, and oil-water viscosity ratio. Reducing oil-water contact angle and oil-water viscosity ratio and increasing oil-water interfacial tension are conducive to the imbibition process. The increase in injection water temperature is usually beneﬁcial to the occurrence of the imbibition. Moreover, the actual core structure imbibition degree is often lower than that of the ideal core structure. The inhomogeneous distribution of rock particles has a signiﬁcant inﬂuence on imbibition. This study provides microscale theoretical support for seeking reasonable injection velocity, pressure gradient, injection temperature, and well-shutting time in the ﬁeld process. It provides a reference for the formulation of ﬁeld process parameters.


Introduction
A large number of artificial fractures are often produced by large-scale hydraulic fracturing in the process of tight reservoir exploitation. Moreover, the reservoir may develop various micron-scale natural microfractures. In the process of water injection to supplement reservoir energy such as water huff-n-puff, oil-water mass transfer and exchange occur between matrix and complex fractures, including oil-water imbibition. Due to the tight reservoir's small pore throat radius, the capillary force is strong, and the imbibition oil recovery has become a critical oil recovery mechanism [1]. The temperature of injection water is often different from that of the reservoir. The temperature directly impacts the fluid properties and the interaction of rock and fluid, thereby affecting the imbibition process. Fluid flow drives heat transfer, and temperature changes affect the properties of the fluid. Therefore, it is necessary to study the coupling of mass and heat transfer of multiphase fluid in porous media at the microscale.
The key is to determine the moving interface between different phases in microscale multiphase heat transfer simulation in porous media. Scholars have simulated it in different ways. The microscale physical experiment simulation [2][3][4][5][6] first received attention.
However, since the limitations of the microscopic physical experiment process, underground reservoir environments (such as high temperature and high pressure) cannot be fully simulated. In recent years, with micro-seepage theory and computer performance improvement, the numerical simulation method based on micromodel has also developed rapidly. The pore network model, lattice Boltzmann method, and interface capture method are widely used in numerical simulation of microscopic rock seepage.
The pore network model [7][8][9] usually uses a simple geometry pore network, which is insufficient to solve the complex pore structure fluid flow. The lattice Boltzmann method [10][11][12][13][14] uses the discrete medium model to simulate fluid motion. Interface capture methods based on Navier-Stokes equations include the volume-of-fluid method (VOF), the level-set method (LSM), and the phase-field method (PFM). The method treats fluid as an uninterrupted whole. It calculates the parameters of fluid motion by calculus, which can deal with the two-phase interface change under complex pore geometry. It is more suitable for the simulation of micro two-phase flow in the complex pore structure. The volume-of-fluid method [15,16] can ensure mass conservation and determine the interface location by solving the volume fraction function. However, the treatment effect is insufficient when the physical quantity of the phase interface changes suddenly. The level-set method [17][18][19] is essentially a mathematical method that uses the smooth function to approximate the motion at the interface and has high calculation accuracy but cannot guarantee quality conservation. In the phase-field method [20][21][22][23][24], based on the physical method, the interface parameters between two phases are combined into the system free energy function. The sharp interface is treated as a diffusion interface layer with a certain thickness. Therefore, the phase-field method can capture more actual physical phenomena in a more reasonable calculation time than the level-set method, which solves the interface problem in multiphase flow and ensures the system's total energy minimization. Amiri [25] compared the level-set method and the phase-field method. It is considered that the phasefield method is more realistic in dealing with the pressure gradient and fluid profile in complex porous media, and the calculation time is significantly reduced compared with the level-set method.
Scholars have studied the mechanism of oil-water imbibition by different methods at the micro-level. Cai et al. [26] presented a state-of-the-art review of the Lucas-Washburn equation-based models and numerical simulation techniques in capillary-driven flow in porous systems. Rangel-German et al. [2] studied the relative flow of fluid between the matrix and the fracture through the glass etching model. The experiment found that when the injection rate was slow, water entered the matrix through imbibition and discharged oil into the fracture. Akbarabadi et al. [27], Mitchell et al. [28], and NMR (Nuclear Magnetic Resonance) were used to study the pore structure of tight core on spontaneous imbibition. The results show that the reservoir maximum connected to pore throat radius is positively correlated with imbibition recovery, and the specific surface is negatively correlated with imbibition recovery. Wang et al. [29] established a pore network model to study the two-phase permeability problem and obtained the phase permeability curves under different contact angles and capillary numbers. Hatiboglu and Babadagli [6] used the lattice Boltzmann method to study the spontaneous imbibition at the pore scale and found that this method has limitations for simulating non-miscible flow with a high viscosity ratio. Zhu et al. [30] studied the change of dynamic wetting angle in a single capillary by the phase-field method. Rokhforouz et al. [31] studied the mechanism of reverse fluid imbibition in heterogeneous porous media by the phase-field method, and regular circles represent rock particles. The model is large and does not consider the influence of temperature change. Jing et al. [32] compared the glass etching experiment with the phase-field method and concluded that the two methods' recovery efficiency and residual morphology were the same. Furthermore, they analyzed the mechanism of surfactant that can enhance recovery. Amiri et al. [21] studied the non-isothermal water flooding process in two-dimensional porous media using a conceptual model and analyzed the temperature distribution in porous media after hot water injection and did not study its influence temperature on the distribution of oil-water two-phase.
Spontaneous imbibition is an important oil recovery mechanism in fractured waterwet reservoirs [33,34]. In this paper, the theoretical and actual models of non-isothermal oil-water imbibition flow in the fracture-matrix of tight reservoirs are established and solved using the finite element method based on the phase-field method. The theoretical heterogeneous model was used to study the mechanism of temperature affecting the imbibition process by changing the oil-water viscosity ratio, contact angle, and interfacial tension. The actual rock model was used to study the oil-water distribution and imbibition effect at different injection temperatures. The effect of temperature change on the imbibition process of tight oil was discussed from a microscopic scale.

Model Geometry and Assumptions
There are many artificial fractures and microfractures in the reservoir after largescale fracturing. In the well-shutting stage, most of the fractures are filled with water or fracturing fluid. Injecting a lot of low-temperature liquid makes the matrix temperature change around the fracture. In this process, imbibition will occur under the guidance of capillary pressure, and the effect of temperature change in the whole process cannot be ignored. In this paper, we first extract the two-dimensional matrix-fracture system, as shown in Figure 1. In order to facilitate the establishment of the model, the following assumptions are made:
Local thermal equilibrium assumption between fluid and rock; 3.
The seepage process is two-dimensional horizontal flow, ignoring the influence of gravity; 4.
Changes in fracture width caused by mechanical factors are not considered.
We established two models to study the problems proposed in this paper: Model I: The first model is an ideal one-sided reverse imbibition heterogeneous porous medium model, as shown in Figure 2. The rectangular area above the model represents the fracture, and the initial fracture is filled with water phase. The matrix area below the fracture is filled with the oil phase. The triangular element is used to mesh the model. Model II: We selected a scanning electron microscope image of a dense matrix in reference [35] to simulate the imbibition process under an actual rock structure ( Figure 3). The morphological operation was used to eliminate fine particles, separate fine-pore throat, and smooth matrix boundary. We used MATLAB software (MathWorks, America) to write code to binarize the image, remove noise, and finally obtain the rock particle real distribution. In the actual model, the shape of rock particles is irregular and irregularly distributed, and a rectangular area is artificially constructed above the original configuration to represent fracture. The simulation parameters of the above two models are unified in Table 1. The equivalent permeability of single-phase water passing through the physical model is calculated to represent the model flowability according to Equation (1). We believe that there is still a small pressure gradient in the fracture during the well-shutting stage after fracturing. We assume that the fracture pressure gradient is 0.001-0.01 MPa/m, and the flow in the fracture meets the fracture flow law between the fixed plates. The average flow velocity is calculated by Equation (2) as the boundary condition of the fracture. The simulation results and analysis will be shown in Sections 3 and 4.
where Q is the water phase flow obtained by integrating the water phase flow velocity at the outlet after reaching the steady-state, m 3 /s; µ w is the viscosity of water, Pa·s; H is the width of porous media, m; and ∇p is the pressure gradient, Pa/m.

Governing Equations
The phase-field method can capture the phenomena related to viscous force and capillary force in a reasonable calculation time. It solves the interface problem in multiphase flow. It ensures that the system's total energy (the combination of surface energy and fluid energy) is minimized, so it has specific physical significance. The phase-field method's main idea is to treat the transition interface of two-phase motion as a thin layer with a width, not zero. A phase-field constant ϑ is introduced to characterize different phases and phase interfaces without directly tracking the two fluids interface changes. ϑ = 1 represents the oil phase, ϑ = −1 represents the water phase, and −1 < ϑ < 1 represents the oil-water interface. Simultaneously, the fluid properties corresponding to each phase-field constant can be expressed as follow: where n can represent the physical properties of oil-water two-phase, such as viscosity, density, specific heat capacity, thermal conductivity, and others. Thus, the smooth transition of fluid properties at the interface is actualized. The control equation mainly includes the Cahn-Hilliard equation controlled by phasefield variables, the Navier-Stokes equation describing mass flow, and the energy conservation equation describing heat transfer. The Cahn-Hilliard convection-diffusion equation governing the change of interface shape is introduced by constructing mixed free energy of interface with constant phase-field [36,37].
where u represents the fluid velocity, m/s; γ represents the mobility; γ c represents the characteristic mobility; ε represents the capillary thickness to characterize the interface thickness, m; G represents the chemical potential of the system and λ represents the mixed energy density; and σ represents the oil-water interfacial tension, N/m.
The source term of Equation (4) is a fourth-order term. ψ is introduced to decompose the fourth-order partial differential equation into two second-order partial differential equations to facilitate calculation [38].
The chemical potential of the system can be rewritten into: Equation (3) can be rewritten to Two critical parameters affecting the accuracy of interface simulation are mobility γ and interface thickness ε. For mobility γ, on the one hand, it must be large enough to keep the interface thickness constant to ensure that the phase-field equation is stable; on the other hand, it must be small enough to maintain the flow accuracy. The reasonable value of mobility is proportional to the interface velocity and inversely proportional to the surface tension coefficient. In the simulation, we take γ = 1, to ensure convergence and obtain more realistic physical results.
Another critical parameter is the interface thickness ε, representing the two-phase interface thickness when the phase-field function changes. The thickness of the interface should not be less than the size of the grid element. When the grid size at the interface is larger than the interface thickness ε, the calculation of interface tension will be wrong, and then the calculation error will be generated. Usually, we set the interface thickness to the size of the maximum grid at the interface. Amiri et al. [25] provided a better convergence at the interface when the mesh size h = 0.8. He also defined the Cahn number (C n = ε/L c ), where L c is the average particle size of porous media. When C n = 0.03, the phase-field model has good convergence and the level-set method needs C n = 0.01, So the level-set method needs a more fine interface grid size.
The surface tension is equivalent to the product of the gradient of the phase-field constant and the chemical potential, and it is introduced into the Navier-Stokes equation as a volume force.
where p is the pressure, Pa; u is the fluid velocity field, m 3 /s; ρ is the density, kg/m 3 ; and µ is the viscosity, Pa·s. The energy conservation equation is based on the local thermal equilibrium assumption, which ignores the convective heat transfer between the surface of rock particles and the fluid. It considers that the temperature of the two parts is equal.
where T is the temperature,°C; C p is specific heat capacity, J/(kg·K); and k is the thermal conductivity, W/(m·K).
The physical quantities ϑ, ψ, u, p, and T are unknowns. For the surface of rock particles, the boundary conditions considering oil-water contact angle are as follows [38,39]: where n represents the unit normal and θ c represents the wetting angle. Further, it is necessary to establish the corresponding relationship between temperature and relevant parameters to study the effect of temperature on imbibition. This paper considers three mechanisms: the relationship between temperature and oil viscosity [µ o (T)], oil-water wetting angle [θ c (T)], and oil-water interfacial tension [σ(T)].
Define dimensionless time as follows: where u inj is the injection rate, t is the time, and L f is the length of the fracture. Defined dimensionless capillary number [21]: COMSOL Multiphysics is a powerful multi-physical field simulation software. We use this software to mesh and solve the finite element problem. The backward difference method is used to discrete the time, and the software's PARDISO direct solver is used to solve the problem. This method's advantage is that the memory is small, and the specific instructions can refer to the software manual [40].

Phase-Field Method Verification
This section verifies the phase-field method applicability in oil-water two-phase flow at the tight oil micro-nano pore scale. The flow in the 10 µm × 1 µm capillary is simulated. The capillary channels are filled with oil in the initial state. Then the water phase enters from the left side. As shown in Figure 4, the blue region represents the water phase, and the red region represents the oil phase. Only capillary force and viscous force are considered. The simulation results are compared with the theoretical formula proposed by Washburn and Lucas [41] (Equation (17)).
where µ w and µ o are the viscosity of water and oil, Pa·s; x 0 is the initial curved surface position, m; L is the capillary length m; r is the liquid surface radius, m; σ is the interfacial tension, N/m; θ is the wetting angle; and t is the time, s. As shown in Figure 5, when oil and water viscosity are 0.005 Pa·s and 0.001 Pa·s, the interfacial tension is 0.01 N/m, and the wetting angle is 30 • , 45 • , and 60 • , the simulation results are in good agreement with the calculation results of the theoretical formula. It is verified that the phase-field method can accurately simulate the two-phase flow process in the micro-nano pore scale.

Thermal Simulation Verification
We assume that the local thermal equilibrium between the matrix and the fluid is satisfied and establish a single fracture heat transfer model, as shown in Figure 6. The cold water is injected with a constant velocity at one end of the fracture, and the analytical solution is given in Appendix A. Model parameters are in Table 2.  Moreover, we used COMSOL Multiphysics to simulate this process. The temperature distribution after one day (86,400 s) of injection is shown in Figure 7. It can be seen that the calculation results of COMSOL Multiphysics are consistent with the analytical solution at different grid sizes, which also shows the accuracy of this research method. As can be seen from Figure 8, the smaller the grid size is, the closer the COMSOL Multiphysics simulation results are to the analytical solution.

Analysis of Tight Oil Matrix-Fracture Imbibition Mechanism
The imbibition can be divided into forwarding imbibition and reverse imbibition. The flow direction of oil and water in forwarding imbibition is the same, and the flow direction of reverse imbibition is opposite to that of oil. In fractured reservoirs, forwarding imbibition mainly occurs in microfractures, which is the process of water entering different scales of fractures and displacing oil. Reverse imbibition mainly occurs in the tight matrix, manifested as water entering small pores first and discharging oil through large pores. We show these two mechanisms in Figure 9. The left side is the forwarding imbibition mechanism, and the right side is the reverse imbibition mechanism. The change from t 1 to t 4 is spontaneous, and the capillary force is the main driving force. We calculate the dimensionless number of criteria to describe physical variables in microscale flow. The Eotvos number (Equation (18)) represents the ratio of gravity to surface tension. The Weber number (Equation (19)) represents the ratio of inertial force to surface tension.
D is the characteristic length. Since the order of magnitude of D in tight reservoirs is 10 −6 m, the calculated Eotvos number and Weber number are very small, so we believe gravity and inertial force can be ignored relative to surface tension. The correctness of this assumption is also shown in the simulation example ( Figure 10).

Micro-Seepage Process of Tight Oil Matrix-Fracture
We take the ideal heterogeneous model in Section 2.1 as an example to analyze the imbibition process. The wetting angle is 22.5 • . The oil viscosity is 3 mPa·s. The oil-water interfacial tension is 25 mN/m. The distribution of the oil phase and water phase at different times is recorded.
As shown in Figure 11 as the above analysis of the mechanism, we can see that the flow channel size difference causes the capillary pressure difference. The capillary pressure at the small pore throat is more significant than that at the large pore throat, and then the pressure difference is formed at different outlet ends of the pore. The water phase enters along the small pore throat, and the oil phase discharges along the large pore throat so that the oil phase is replaced by the fracture in the state of small oil droplets. The small oil droplets in the fracture continuously gather into large oil droplets. Then they flow forward along the fracture wall under the action of viscous force (caused by the flow velocity in the fracture) until no small oil droplets are replaced and gradually tend to be in the state of osmotic equilibrium. This is consistent with the phenomenon observed in the experiment [2,6]. As shown in Figure 12, we can divide the imbibition process into the fast imbibition stage (I), slow imbibition stage (II), and imbibition equilibrium stage (III), which provides micro-mechanism support for the optimal well-shut time in the water huff-n-puff process. The injection temperature is set as 20 • C. It can be seen from Figure 13 that the lowtemperature part can spread to the matrix in a short time. The temperature change will affect the viscosity, wettability, and interfacial tension of the fluid, affecting the imbibition process analyzed in Section 4, which provides micro-mechanism support for the optimal injection temperature. In the pressure field distribution and velocity field distribution (Figure 14), it can be seen that small oil droplets have more significant internal pressure, which means that large oil droplets are easier to deform and be displaced. We found that the generated small oil droplets hindered the original water flow. It can be predicted that when the generated oil droplets are fast enough and large enough, a greater viscous force is needed to promote the migration of these oil droplets in the fracture. However, a greater viscous force is not conducive to the imbibition process, and a reasonable pressure gradient or injection rate needs to be determined. We can provide micro-mechanism support for the need to seek reasonable injection temperature, injection rate, and well-shut time in the water huff-n-puff process through the above analysis of the imbibition process.

Results and Discussion
We consider that temperature affects the imbibition process by changing fluid viscosity, oil-water interfacial tension, and oil-water wettability. Firstly, the influence of the wetting angle, the interfacial tension, and the oil-water viscosity ratio on the imbibition process are analyzed using a theoretical heterogeneous model. The purpose of the single-factor analysis is to study the influence mechanism of temperature change on the matrix-fractured oil-water system's imbibition process and the influence trend of each mechanism on the imbibition effect. Furthermore, the influence of temperature on the imbibition process was studied by considering the actual rock structure.

Wetting Angle
Oil-water imbibition only occurs in water-wet reservoirs. Wettability is an important factor affecting the imbibition process. The wettability of rock changes at different temperatures. Different fluid-rock systems often show different laws for the effect of temperature on the wetting angle. Hamouda [42] experimentally measured that with the increase in temperature, the wetting angle decreased, and calcite changed from oil-wet to water-wet. The higher the temperature was, the stronger the rock hydrophilicity was. Wang [43] improved the suspended droplet interfacial tension meter. The sandstone system's contact angle increased with the increase in temperature, and the contact angle of the carbonate system decreased with the increase in temperature. The experimental results of Hjelmeland et al. [44] show that the wettability of rock changes from oil to water with increasing temperature. Rajayi et al. [45] measured through experiments that for quartz/water/bitumen systems, the rock will change from moisture to neutral moisture with the increase in temperature. H. Jabbari et al. [46] believe that hot water injection can change the matrix from the neutral or oil-wet state to the water-wet state, thereby enhancing oil recovery.
Therefore, we believe that the influence of temperature on wettability is different under different rock-fluid systems. In this part, only the mechanical analysis is performed, and the wetting angle is used to represent the wettability of rock. Setting static wetting conditions wetting angles are 90 • , 60 • , 45 • , 30 • , and 22.5 • represent neutral wet, weak water-wet, water-wet, strongly water-wet, and extremely water-wet conditions.
We compared the oil-water distribution after imbibition equilibrium. It can be seen from Figure 15 that under neutral wetting, the oil-water interface only moved slightly to the fracture. The recovery degree of the strongly water-wet reservoir is higher than that of the weak water-wet reservoir. With the decrease of wetting angle from 90 • to 30 • , the area of water displacement oil expands. Furthermore, although the overall suction area becomes more extensive with the decrease of wetting angle from 30 • to 22.5 • , more oil remains in the middle (the area marked by the pink circle in Figure 14). It can be inferred that the uneven distribution of rock particles is also an essential factor affecting the imbibition process. In the process of imbibition, the water phase mainly advances along the rock wall. The more water-wet the rock particles are, the greater the effect of uneven distribution on imbibition is. It can also be seen from Figure 16 that the more water-wet the rock is, the longer time it takes for the imbibition to reach a stable state, and the greater the degree of imbibition recovery of the matrix is.

Oil-Water Interfacial Tension
Hjelmeland et al. [44] obtained through experiments that the glass cell oil-water interfacial tension decreased with temperature increasing; in contrast, it increased in the highpressure cell. Sajjad Esmaeili et al. [47] also measured poly-alpha-olefin (PAO) oil-water interfacial tension decreased significantly with increasing temperature (23 • C correspondings to 41.1 mN/m, 185.2 • C correspondings to 20.9 mN/m). This section only explores the mechanism, setting the oil-water interfacial tension values are 5 mN/m, 25 mN/m, and 50 mN/m. The corresponding dimensionless capillary numbers are 2 × 10 −4 , 4 × 10 −5 , and 2 × 10 −5 . Figure 17 is the oil-water distribution at t D = 8.6538. It can be seen from Figure 18 that the greater the oil-water interfacial tension is, the greater the imbibition rate and degree are. This is because increasing the interfacial tension will lead to increased capillary force, which is the imbibition process's main driving force.   We compared the oil droplets under different interfacial tensions at the same position ( Figure 19). It can be seen that with the decrease in interfacial tension, oil droplets become more prone to deformation, and it is easy to generate some smaller oil droplets into the fracture. When t D = 8.6538, Ca = 4 × 10 −5 and Ca = 2 × 10 −5 entered the imbibition equilibrium stage, and the oil-water interface no longer changes. However, in the case of Ca = 2 × 10 −4 , the imbibition process is still ongoing. We recorded the oil-water distribution at different times in the Ca = 2 × 10 −4 ( Figure 20). About 56% of the oil in the matrix can be extracted when t D = 48.0769. The imbibition rate reduces with the interfacial tension reduction, resulting in a long time to achieve the imbibition balance and produce more heterogeneous residual oil.

Oil-Water Viscosity Ratio
The oil-water distribution at t D = 8.6538 is shown in Figure 21. With the decrease in crude oil viscosity, the spread range of water is more reflected on the left side of the model. For example, when µ 0 /µ w = 3, a fingering phenomenon is on the left side, and the oil-water distribution on the right side changes little.  Figure 22 shows that the oil-water viscosity ratio effect is not evident in the early imbibition stage. The imbibition rate decreases with the decrease in the oil-water viscosity ratio because the viscous seepage resistance increases with oil viscosity.

Simulation Results and Analysis of the Actual Example
This section studies the imbibition process of an actual rock structure, considering the influence of different temperatures on the oil viscosity, oil-water interfacial tension, oil-water contact angle, and the specific fitting expression (Equations (20)- (22)). Figure 23 shows the variation of different parameters with temperature. The data points come from the results of physical experiments or numerical simulation experiments. The injection temperature of 120 • C is taken as an example to illustrate the change of temperature field. Figure 24 shows the temperature change of the fluid region and the rock region. It can be seen that the area affected by the injection temperature is increasing until the whole matrix is affected. It can also be seen from the change process of the temperature field that the thermal convection between fractures and matrix caused by seepage is a slow process relative to thermal conduction. Because the given oil viscosity, oil-water interfacial tension, and oil-water contact angle are only affected by temperature, each grid's fluid properties and fluid-rock interface properties change with time. For example, the grid on the surface of a rock in the region is water-wet at 100 • C, and with the increase in temperature, the surface of the rock becomes more water-wet, so it is more prone to seepage. On the contrary, when the temperature decreases, the rock surface gradually tends to be neutral wetting, resulting in this area not continuing to imbibition. Figure 25 shows the average temperature of the model at different times. At t D = 4.9652, the average temperature of the model is close to the injection temperature.  The oil-water distribution in the whole process is analyzed. It can be seen from Figure 26 that when the injection temperature is 60 • C, with the gradual advancement of the low-temperature zone, the imbibition process tends to stop, and only a small part of the oil phase enters the fracture and flows with water. When the injection temperature is 120 • C, the oil droplets entering the fracture become more easily deformed, caused by decreased interfacial tension caused by increased temperature. As the rock becomes more water-wet, the degree of imbibition becomes greater. When the injection temperature is 180 • C, more oil phase is replaced. Finally, a large oil droplet is formed in the fracture. It can be seen from Figure 27 that when the parameters given in this paper change with temperature, the higher the injection temperature is, the more conducive it is to imbibition. We also found that under the actual complex pore throat connection, the imbibition phenomenon was not as apparent as the ideal model, which was similar to the conclusion obtained by Farzaneh Ramazari [48] using the molecular experiment method.

Conclusions
In this work, the Cahn-Hilliard phase-field, Navier-Stokes equation, and energy conservation equation are coupled to establish the ideal inhomogeneous rock pore structure model and the actual rock pore structure model. The finite element method is used to solve the model. The oil-water imbibition mechanism of the microscale matrix-fracture system under the influence of temperature is studied. The main conclusions are as follows: 1.
The oil-water imbibition phenomena in the matrix-fracture system mainly include the reverse imbibition in the matrix and the forward imbibition in the micro-fracture. The capillary force is the main driving force. At the micro-nano scale, the influence of gravity and the imbibition process can be ignored. The seepage of the actual pore throat model is weaker than the ideal pore throat model.

2.
The oil-water imbibition process of matrix-fracture at the microscale can be expressed as the fracture is first filled with water. Under the action of capillary force, the water phase enters the matrix from small pores. Due to the different pore throat sizes, the pressure difference is formed. Under pressure difference, the oil phase is displaced to the large pores and enters the fracture to form small oil droplets. The small oil droplets gradually gather into large oil droplets, and positive imbibition occurs in the fracture. The oil phase flows along the fracture with the water phase to complete the matrix-fracture oil-water imbibition replacement.

3.
Wettability, oil-water interfacial tension, oil-water viscosity ratio, and inhomogeneous distribution of rock particles affect the matrix-fracture imbibition process. Reducing oil-water contact angle and oil-water viscosity ratio and increasing oil-water interfacial tension is beneficial to the imbibition process.

4.
The heat conduction process is the primary heat transfer mechanism at the microscale. When the injection temperature is different from the formation temperature, it will soon spread to the imbibition area. Therefore, the influence of temperature change in the imbibition process cannot be ignored. Under the conditions given in this paper, the increase in injection temperature is conducive to imbibition.
It is necessary to find a reasonable injection rate, pressure gradient, injection temperature, and well-shutting time conducive to the field process imbibition process. In this regard, this study provides theoretical support at the microscale. It is proved that the phase-field method is reasonable in the study of oil-water two-phase microscale flow. In future work, quantitative research and characterization of influencing factors need to be further studied and solved.