Shale oil is generally stored in organic-rich shale with permeability of less than 0.1×10^-3 μm² and porosity of less than 10% ^[1-2]. One of the main methods to improve shale oil recovery is to greatly increase the drainage area by adopting long horizontal wells and multi-stage and multi-cluster volume fracturing. Field production data indicate that even with the adoption of volume fracturing and long horizontal well, a large number of shale wells are still not commercially viable ^[3]. Fracturing-soaking technology provides a new idea for shale oil development^[4]. The Changqing, Qinghai, Daqing, Jilin and other oil fields have adopted the production management system of "fracturing, soaking and re-production" in some shale oil blocks successively, and achieved good stimulation effect. At present, most shale oil wells in China are shut in for 14-60 d after fracturing, and the determination of shut-in time mostly, without scientific and reasonable basis, depends largely on the experience of operators ^[5-6]. Establishing the optimization model and method of shut-in time after fracturing is the key to the implementation of the fracturing and soaking, and the theoretical and engineering issue needs to be solved for further improving shale oil recovery.

The purpose of fracturing is to reconstruct the reservoir to obtain larger fracture distribution area and increase the high-permeability volume space for oil and gas flow. The purpose of shut-in is to increase the flow energy of formation fluid and enhance the imbibition effect of oil-water displacement by utilizing a large amount of fracturing fluid retained in the reservoir based on the coupling effect of fluid flow, imbibition principle and chemical reaction, and finally achieve the purpose of increasing crude oil production ^[7]. Although simple in operation, shut-in is very complex in core physical mechanism. Lee et al. ^{[8⇓⇓⇓⇓-13]} mainly studied the imbibition behavior dominated by gravity and capillary pressure at the core scale and field scale by lab experiment and numerical simulation. They concluded that the change of rock wettability and complexity of reservoir reconstruction were important factors affecting the stimulation effect of soaking. Li et al. ^{[14⇓⇓⇓⇓-19]} considered the laminar flow caused by chemical potential and the influence of osmotic pressure besides the capillary pressure dominated imbibition on the stimulation effect of soaking. They found that the difference of salinity inside and outside the matrix in the process of soaking was the fundamental reason for generating osmotic pressure and promoting crude oil displacement. Bui et al. ^[20-22] also mentioned in their study that the elastic energy of the reservoir would be supplemented due to the filtration of fracturing fluid into the reservoir, which might also lead to the increase of productivity after shut-in. If the reservoir is deemed oil wet, adding surfactant to the fracturing fluid can further improve the wettability of the reservoir through the fluid flow during shut-in to increase production. Generally speaking, most of the researchers studied the mechanisms of imbibition oil recovery during shut-in, and rarely examined the optimal shut-in time after fracturing. Although carrying out a detailed study on the optimal soaking time, Wang et al.^[23] didn’t consider the influence of osmotic pressure on imbibition in the model and the multiphase flow in hydraulic fractures, and the physical model they established wasn’t a complete filed scale model.

This work aims to establish an integrated shut-in time optimization model of fracturing, soaking and production for shale oil reservoir, and puts forward a soaking time optimization method with the goal of maximizing productivity considering the effects of capillary pressure, chemical potential, osmotic pressure and hydraulic fractures. The model has been verified by an application example and commercial software. In the integrated simulation process of fracturing-soaking-production, the key production-increasing mechanism of soaking has been analyzed from multiple angles, and the influence laws of different factors on the optimal soaking time of shale oil development have been investigated emphatically, so as to guide the decision-making of soaking time after fracturing of shale oil horizontal wells and provide a theoretical basis for further optimization of shut-in system.

Imbibition production is an important mechanism for enhancing oil recovery of shale oil reservoir. Besides capillary pressure and gravity, chemical osmotic pressure, membrane properties of matrix wall and fluid exchange area between matrix and fracture are also important factors affecting imbibition during shut-in. Regardless of fluid damage, at different shut-in times, the displaced oil amount by imbibition and cumulative production are different. With the increase of shut-in time, the cumulative production or recovery increase to a certain extent.

Based on the established fracturing-soaking-production integrated model, a shut-in time optimization method aiming at maximizing productivity and recovering cost fastest has been proposed in this paper. The specific operation steps are as follows.

(1) According to the established model, the recovery performance under a certain set of engineering parameters, geological parameters and a certain shut-in time is obtained.

(2) Change the shut-in time and repeat step (1). Under the same geological parameters and engineering parameters, the relationship between cumulative output and production time under different shut-in time can be obtained, as shown in Fig. 5a.

(3) Based on the results of step (2), sort out the differences in cumulative production between different shut-in time schemes and non-soaking scheme at a certain production time (generally select the time point with large difference in cumulative productions), so as to obtain the relationship between the cumulative production increment and shut-in time shown in Fig. 5b.

(4) Based on Fig. 5b, the optimal shut-in time is picked out. In order to select the optimal shut-in time quantitatively and ensure that the cumulative production increment basically reaches the peak value within the shortest shut-in time, the first derivative of the relation curve between the cumulative production increment and shut-in time can be calculated. When its value is 0.001 t/d, that is, the slope of a point on the curve is 0.001 t/d, the corresponding shut-in time of this point is considered as the optimal shut-in time. As shown in Fig. 5b, the slope at point A is 0.001 t/d, and the corresponding optimal shut-in time is 30 d. At this point, the flow of oil, water and solute in the system is basically in equilibrium state, and the cumulative production increment is reaching the peak. If the well remains shut-in for longer, the increase of cumulative production increment is small. The shut-in time of 30 d can achieve the goal of the fastest cost recovery.

The optimal shut-in time is the time for the multi- phase fluid flow and energy distribution to re-balance again in reservoir after fracturing and is a key design parameter in the soaking decision. According to the integrated fracturing-soaking-production model and shut- in time optimization method presented in this work, combined with the fracturing pressure data of shale oil well C in Jimsar sag, a single-stage with 5-fracture-cluster physical model was established, as shown in Fig. 6. The 1-5 clusters of fractures in the model were 142, 82, 120, 82 and 142 m in half-length respectively, and each cluster was 40×10^-3 μm²∙m in average conductivity. The three continuous processes, fracturing, soaking and production, were simulated continuously. First, a short time injection at high injection rate was conducted to simulate the fracturing process, then the source term was set at zero to start well soaking, and finally the well was opened for production. The main input parameters in the simulation are shown in Table 5, and the relative permeability curve and capillary pressure curve are shown in Fig. 3.

Figs. 7 and 8 show the distribution of pore pressure and oil saturation at different moments in fracturing, soaking and production stages respectively. It can be seen that fracturing is a physical process similar to transient water flooding. With the increase of fracturing time, oil saturation near fractures gradually decreases while the pore pressure gradually increases. In the soaking stage, influenced by capillary pressure, chemical potential and osmotic pressure, the oil saturation near the fractures gradually recovers, and a large amount of crude oil in the matrix is displaced into the fractures, so the oil production is high at the beginning of the well opening after soaking. With the increase of shut-in time, pore pressure tends to a certain equilibrium value. Under closed boundaries, the pore pressure when the reservoir reaches equilibrium is greater than the initial pore pressure. In production stage, the simulation results of pore pressure and oil saturation show the same trends as in conventional production process.

By changing the shut-in time, the production performance curves at different shut-in moments were simulated (Fig. 9a). Taking the cumulative production of 300 d into production without shut-in as the reference, the relationship between the soaking time and the cumulative production increment is shown in Fig. 9b. Fig. 9a shows that the initial daily production and cumulative production increased significantly after shut-in, but the increase of the cumulative production decreased with the increase of production time. Fig. 9b shows that with the increase of shut-in time, the cumulative production increment increased first and then gradually approached a stable value. In order to reasonably control shut-in time, produce crude oil and recover the cost as soon as possible, according to the shut-in time optimization method, the point with curve slope of 0.001 t/d in Fig. 9b was found, and the corresponding value of 50 d was the optimal shut-in time.

Based on the integrated fracturing-soaking-production model and the physical model shown in Fig. 6b, the influences of matrix permeability, porosity, membrane efficiency, capillary pressure multiple, displacement rate, total injected liquid volume and fracture length on the optimal shut-in time were studied. When analyzing the influence of a parameter, only different values were set for this parameter. The values of other parameters were the same as those in the above example.

With the increase of shut-in time, both the initial production and cumulative production increment of shale wells increase rapidly at first and then tend to a stable value. The shut-in time corresponding to the inflection point in the trend of cumulative production increment is the optimal shut-in time.

The optimal shut-in time has a nonlinear negative correlation with matrix permeability, porosity, capillary pressure multiple and fracture length, a nonlinear positive correlation with membrane efficiency and the total injected liquid volume, and a nearly linear positive correlation with displacement rate. The factors in descending order of influence on the optimal soaking time are the total injected liquid volume, capillary pressure multiple, matrix permeability, porosity, membrane efficiency, fracturing fluid salinity and displacement rate.

c_sm, c_sf—solute mass fraction in matrix and fracture, %;

c_sf,in,pro—solute mass fraction of sink term in fracture, %;

c_swmf—solute mass fraction in the unit where fracture and matrix intersect, %;

Δc_smf—difference of solute mass fraction between matrix and fracture, %;

C_m, C_f—compressibility coefficients of matrix and fracture porosity, Pa^-1;

D_eff—effective volume diffusion coefficient in porous media, m²/s;

D_m—altitude, m;

E_m, E_f—stress sensitivity coefficients of matrix and fracture permeability, Pa^-1;

F—level of significant difference;

K_m, K_f—permeability of matrix and fracture, m²;

K_0,m, K_0,f—initial permeability of matrix and fracture, m²;

n—fracturing stages;

p₀—initial pore pressure, Pa;

p_c—capillary pressure, Pa;

p_of, p_wf—pressures of oil phase and water phase in fracture, Pa;

p_om, p_wm—pressures of oil phase and water phase in matrix, Pa;

Δp_omf, Δp_wmf—pressure differences of oil phase and water phase between matrix and fracture, Pa;

Δp_op—osmotic pressure difference, Pa;

P—detection level;

q_of, q_wf—source and sink terms of oil and water phases in fractures, s^-1;

q_omf, q_wmf, q_smf—exchange terms of oil phase, water phase and solute between matrix and fracture, s^-1;

R—gas constant, J/(mol·K);

S_of, S_wf—saturation of oil phase and water phase in fracture, %;

S_om, S_wm—saturation of oil and water phases in the matrix, %;

t—time, s;

T_m—temperature of matrix, K;

T_omf, T_wmf—conductivity coefficients of oil phase and water phase between matrix and fracture based on embedded discrete fracture model, m³;

V_m—grid cell volume, m³;

V_w—partial molar volume of water, m³/mol;

X, Y—rectangular coordinate system, m;

γ_o, γ_w—API gravity of oil phase and water phase, N/m³;

ε_vm—volume strain, dimensionless;

λ_of, λ_wf—relative fluidity of oil phase and water phase in fracture, (Pa•s)^-1;

λ_om, λ_wm—relative fluidity of oil phase and water phase in matrix, (Pa·s)^-1;

λ_omf, λ_wmf—relative fluidity of oil phase and water phase corresponding to matrix and fracture intersecting unit, (Pa·s)^-1;

ρ_w, ρ_s—density of the water phase and solute, kg/m³;

ϕ_m, ϕ_f—porosity of matrix and fracture, %;

ϕ_0,m, ϕ_0,f—initial porosity of matrix and fracture, %;

ω—membrane efficiency, dimensionless.