The large-scale hydraulic fracturing technology of horizontal wells enables effective development of shale oil. Horizontal wells in shale oil reservoirs are often shut in at first after fracturing, and then are put into production, and field practices have also proved that shut-in after fracturing has certain stimulation effect. During the shut-in process, the fracturing fluid in the hydraulic fracture imbibes into the reservoir matrix and replaces oil and gas, which has become an important stimulation mechanism of shut-in in shale reservoirs^[1]. At present, the shut-in period of shale oil fields in China is mostly 14-60 d. Downhole pressure is not measured in the oil wells during shut-in period generally, but only the wellhead tubing pressure data is recorded. These data has not been well utilized due to low measuring accuracy.

In terms of the evaluation of imbibition and oil replacement effect, many scholars have done physical simulation experiments on the core scale^{[2,3,4,5,6,7,8,9,10]}. The results show that capillary pressure is the main driving force of imbibition, and the complexity of hydraulic fractures and wettability of reservoir matrix are important factors affecting the effect of fracturing fluid imbibition and oil replacement. These research results provide a theoretical basis for developing more complex hydraulic fractures and fracturing fluid formulas more conducive to imbibition and oil replacement. In addition, the results of numerical simulations of imbibition and oil replacement on the reservoir scale^{[11,12,13,14,15,16,17,18,19]} show that the imbibition and oil replacement effect of reservoir matrix can be successfully characterized by setting capillary pressure, which can guide the optimization of production regime during shut-in after fracturing and fracturing operation scheme. In addition to the forward physical simulation experiment and numerical simulation research, scholars tried to use the inversion theory of pressure drop well test interpretation to evaluate the fracturing effect of shale reservoirs^{[20,21,22,23,24,25,26,27]}, but there are few inversion studies on fracturing fluid imbibition and oil replacement in the evaluation of fracture morphology and conductivity. There are two reasons: (1) Previous injection pressure drop models are mostly calculated by analytical solutions, which have the limitation of considering fracturing fluid as water and assuming single-phase water flow in fractures and reservoirs, without considering the leak-off of fracturing fluid into matrix pores to replace crude oil and the consequent changes of two-phase flow and pressure field. (2) Shale reservoirs are tight. Whether small-scale fracturing test or pump-off pressure drop test after primary fracturing, the test time generally lasts for only a few hours, and fractures have not closed in such short time, so the effect of oil replacement by fracturing fluid imbibition cannot be evaluated by inversion method.

In this study, by making use of pressure drop data during long period of shut-in, a pressure drop model and analysis method have been worked out to diagnose the dynamics of oil-water displacement between well, fracture and reservoir in the process of shut-in after fracturing and calculate fracture network parameters, matrix imbibition volume and oil replacement volume. This method has been proved valid through example analysis, and can be used to guide the evaluation of fracturing effect of shale oil horizontal wells and provide a theoretical basis for further optimization of shut-in production regime and fracturing fluid formula.

The assumptions are as follows. (1) The fractured shale reservoir consists of a horizontal wellbore (W), main hydraulic fractures (F), secondary fractures (f) and matrix (m), in which wellbore is the medium for fluid exchange between reservoir and surrounding. (2) The fracturing- fluid injection is performed at multiple stages simultaneously, after which the bridge plugs dissolve completely to allow fluid transport through the wellbore. The connection mode of the four media are as follows: W is connected with F, F is connected with both f and m, and f is connected with m. (3) Oil-water two-phase flow is isothermal. (4) F, f and m are compressible. (5) The capillary imbibition of matrix is considered. (6) One main hydraulic fracture with two identical wings distributed symmetrically on the two sides of the horizontal wellbore is created in each fracturing cluster; and the main hydraulic fractures and secondary fractures are vertical fractures with the height the same as the thickness of the reservoir.

The physical model of multi-stage hydraulic fractured shale oil horizontal well is simplified as the combination of W, F, f and m based on the above assumptions, coupled by the continuities of pressure and flow rate at the interfaces. The grid representation of the physical model is shown in Fig. 1. The x, y and z are used to characterize the three directions of the model. The wellbore is the source-sink term. The grid of bottom-hole is located at the center of the main hydraulic fracture. The main hydraulic fracture is the vertical fracture distributed symmetrically at the two sides of the horizontal wellbore, and represented by refined grids with high conductivity. It is characterized by length, width and height. Secondary fractures are vertical fractures which are orthogonal to the main hydraulic fractures. Porosity and permeability are used to characterize the conductivity of secondary fractures. The shape factor^[28] is used to characterize the density of secondary fractures. The matrix is evenly distributed around the main hydraulic fractures and secondary fractures. The storage capacity and flowing ability of the matrix are characterized by porosity and permeability, respectively. The mass transfer process of oil-water two phases during shut-in is shown in Fig. 1. The pumping rate at the surface is zero after the pumping of fracturing-fluids. The fracturing fluid in the wellbore continues to flow into the main hydraulic fractures under pressure difference. The fracturing-fluids in the main fractures flow into the secondary fractures and leak into the matrix. The fracturing-fluids in secondary fractures enter the matrix under the hydraulic pressure difference and capillary pressure. The oil in the matrix is replaced to the secondary fractures and main hydraulic fractures.

Mass balance equation of water phase is:

Mass balance equation of oil phase is:

Fluid exchanges occur between adjacent media, but they cancel each other out in Eqs. (1) and (2). The water fluxes between two adjacent media are given by:

In consideration of the compressibility of fractures and matrix pores, the stress sensitive equations of porosity and permeability are given as:

In consideration of two-phase flow of water and oil in fractures and matrix, the constraint equation of water saturation needs to be added:

In addition, capillary pressure is not considered in the fractures with high permeability, while the capillary pressure in the matrix is given by:

Assuming that the matrix and fractures have the same initial pressure and the same initial water saturation, which is under the original reservoir conditions before development. Constant bottom-hole pressure was set for fracturing-fluids injection. The porosity compressibility coefficients and permeability compressibility coefficients (Eqs. (7) and (8)) of the fractures and matrix were adjusted to control the pumping rate of fracturing-fluids to ensure that the injection volume of the model is the same as that of the actual case during the fracturing treatment. The shut-in pressure drop modelling was initialized by the pressure distribution and water saturation distribution at the end of the injection, which is given by:

Closed boundary is set to be the outer boundary condition.

Finite difference method was used to discrete the equations above. The difference schemes of time and space were processed by front difference method and central difference method, respectively. The implicit method was used to characterize the unknowns. The semi-implicit method was used to linearize the coefficients of nonlinear terms. The upstream weighted method was used to obtain conductivity values. After substituting the initial and boundary conditions, the coefficient matrix was extracted. The Gauss-Seidel method was used to solve the system of linear equations to obtain the values of the unknown variables at the current time step. The above steps were repeated to obtain the values of the unknown variables at the next time step until the set calculation time was reached.

The model was established based on the geological and operation parameters of a typical fractured horizontal well in Jimsar shale oil reservoir of Xinjiang. The reservoir had an original pressure of 38 MPa, and the length, width and height of 1500 m, 560 m and 42 m, respectively. A horizontal well with a horizontal section of 1200 m in length was drilled in the center. The well was fractured in 30 stages with 3 clusters in each stage (a total of 90 main hydraulic fractures). Each main hydraulic fracture was 140 m in half-length and 8 μm²·cm in conductivity. The secondary fractures had a permeability of 0.1×10^-3 μm² and a porosity of 15%. The reservoir matrix had a porosity of 8.6% and a permeability of 0.007×10^-3 μm². The viscosities of water phase and oil phase were set at 3.6 MPa•s and 10 MPa•s, the densities of water phase and oil phase were set at 1000 kg/m³ and 843 kg/m³ respectively. In this model, oil-water relative permeability and capillary pressure of shale matrix were obtained from core experiments, as shown in Fig. 2. The fracturing fluid pumping process was simulated as a water injection process of 2 h with a total fluid volume of 31 500 m³. After the completion of pumping, the pressure drop of shut-in for 120 d was simulated.

During the shut-in process, oil and water phase crossflow exists between wellbore, main hydraulic fractures, secondary fractures and matrix. Fig. 3 shows the simulation results of the after-flow velocity of fracturing fluid from the wellbore to main hydraulic fractures, the crossflow velocity of fracturing fluid from main hydraulic fractures to secondary fractures in the fracture network system, the leak-off velocity of fracturing fluid from main/secondary fractures to the matrix, and the oil replacement velocity from the matrix to main/secondary fractures during the period of the end of pumping to 120 d into shut-in. The simulation results show that: At the 0.01 d after the pumping stopped, the wellbore after-flow phenomenon was obvious, with the after-flow velocity decreased rapidly and then rose slowly until the 5.7 d of shut-in, but still remained extremely low. The fracturing fluid in the main hydraulic fracture would further flow to the secondary fracture and leak into the matrix during the shut-in process. As the shut-in time prolonged, the cross-flow velocity from the main hydraulic fracture to the secondary fracture decreased constantly, and was surpassed by the leak-off velocity of fracturing fluid from the secondary fracture to the matrix at 0.53 d into shut-in. The crude oil in the matrix was replaced into the secondary fractures from the initial moment of shut-in, and the oil replacement velocity of the secondary fracture began to rise and then dropped slowly after 0.1 d. After shut-in for 1.29 d, the oil in secondary fractures began to flow into main hydraulic fractures, and the oil replacement velocity went up gradually till 64.6 d into the shut-in, after that the oil replacement velocity dropped continuously. After shut-in for 1.29 d, oil in the matrix started to be replaced into the main fractures, and the oil replacement velocity rapidly decreased to zero after 5.7 d.

Fig. 4 shows the simulated cumulative water imbibition volume of matrix and cumulative oil replacement volume of main and secondary fractures with shut-in time. During the shut-in of 120 d, a total of 160 05.77 m³ of fracturing fluid flew into the matrix through main hydraulic fractures and secondary fractures, accounting for 50.81% of the total fracturing fluid pumped in. A total of 8430.23 m³ of shale oil was replaced from the matrix by this part of the imbibed fluid, and 541.52 m³ of shale oil flew into main hydraulic fractures, 7888.71 m³ of shale oil flew into secondary fractures. The simulation results show that the secondary fractures are the main body connecting the reservoir matrix, and the main contributor of fracturing fluid imbibition and oil replacement.

Based on the bottom-hole flow pressure during 120 d of shut-in obtained from simulation, the characteristics of bottom-hole pressure drop were analyzed by using the log-log of pressure drop and derivative of pressure drop defined by Bourdet et al.^[29]Fig. 5 shows the log-log curve of pressure drop and pressure drop derivative during shut-in. It can be seen that the bottom-hole pressure drop increased in the initial stage of 0.53 d of shut-in, and then tended stable. The pressure drop derivative curve generally presents W-shape, showing a rising trend at the beginning of shut-in, two concave lows in the middle, and fast decline at the late stage of shut-in (64.6 d), and the pressure drop derivative dropped to zero at 102.5 d into shut-in.

Based on the water imbibition and oil replacement velocities between various levels of media, the pressure drop is divided into eight dominated flow stages (Fig. 5), which are described according to the sequence of shut-in time below.

① Wellbore after-flow-dominated stage (0-0.04 d): The fracturing fluid stored in the wellbore during the pumping stage continued to flow into main fractures, which is shown as coincidence of the curves of pressure drop and pressure drop derivative with the slope of 1.

② Inter-fracture channeling-dominated stage (0.04- 0.53 d): the fracturing fluid channeling existed between main hydraulic fractures and secondary fractures, and the crossflow velocity exceeded the bottom-hole after-flow velocity. The main hydraulic fractures closed, and the slope of pressure drop derivative curve is 1/2.

③ Fracture network leak-off-dominated stage (0.53-1.29 d): the bottom-hole after-flow and inter-fracture crossflow velocities fell fast. The leak-off velocity from fracture network to matrix began to exceed the inter-fracture crossflow velocity, so secondary fractures started to close, the oil replacement velocity from matrix and the fracture network leak-off velocity dropped synchronously, the pressure drop curve rose slowly, and the pressure drop derivative curve showed negative slope.

④ Fracture network storage stage (1.29-5.70 d): The pressure spread to the fracture network boundary and the fracture network pressure began to decline. The oil in the matrix and secondary fractures began to flow into the main hydraulic fractures, and the velocity of oil flowing into main fractures increased, while the leak-off velocity decreased steadily, which is reflected as positive slope of the pressure drop derivative curve.

⑤ Matrix imbibition and oil replacement stage (5.7-8.1 d): Wellbore after-flow velocity reduced to a very low value, the crossflow velocity from main fractures to secondary fractures and the leak-off velocity from main fractures to the matrix of fracturing fluid maintained at low levels, this stage was dominated by oil replacement from the matrix by imbibing fracturing fluid from the secondary fractures, the leak-off velocity of fracturing fluid from secondary fractures fell fast, the oil replacement velocity increased rapidly, the fracturing fluid leak-off velocity was still greater than the oil replacement velocity, which is represented by negative slope of the pressure drop derivative curve.

⑥ Steady increase stage of oil-water exchange efficiency (8.1-64.6 d): The fracturing fluid leak-off velocity and oil replacement velocity of fracture system rose synchronically, and the oil replacement efficiency gradually stabilized at the optimal value, which is represented by positive slope of the pressure drop derivative curve.

⑦ Oil-water replacement efficiency decline stage (64.6-102.5 d): pressure spread to the control boundary, fracturing fluid leak-off velocity and oil replacement velocity dropped synchronously, but the oil replacement velocity dropped faster, resulting in reduction of oil-water displacement efficiency, which is manifested as rapid decrease of pressure drop derivative.

⑧ Oil-water replacement equilibrium stage (102.5- 120.0 d): The bottom-hole pressure stopped dropping, the pressure of the system reached a balance, the oil-water displacement velocity dropped to the lowest, which is represented by zero of pressure drop derivative.

In order to compare the influence of secondary fractures and capillary imbibition on the pressure drop during shut-in period, the porosity and permeability of secondary fractures in the basic model were set the same as those of the matrix to simulate the situation without secondary fractures. The matrix capillary pressure was set as zero to simulate the case without spontaneous matrix imbibition. The pressure drops of the two groups of comparison models above under the same shut-in system were simulated respectively, and the obtained pressure drop characteristic curves were compared with that of the basic model, as shown in Fig. 6. Comparison of the grey pressure drop derivative curve of the case without secondary fractures with the blue curve of basic model shows: if there are only main hydraulic fractures in the fractured reservoir, the leak-off stage ③ and storage stage ④ dominated by secondary fractures would not exist. In addition, when there are no secondary fractures in the fractured reservoir, and the contact area between main hydraulic fractures and the matrix is limited, so the oil-water displacement efficiency would begin to decrease rapidly, which is shown by the short duration of flow stages ⑤ and ⑥, and the earlier coming of stage ⑦. Hence, the pressure drop derivative curve no longer takes the W-shape. Comparison of the yellow pressure drop derivative curve of the case without capillary imbibition with the blue curve of basic mode shows that: if the fracturing reservoir is oil-wet, there will be no capillary imbibition, and the fracturing fluid will flow into the matrix from the fracture wall only under the action of hydraulic pressure difference. In this case, the stage ⑤ and stage ⑥ controlled by matrix imbibition will be smaller in pressure drop derivative values and shorter in duration, consequently, the second concave area of the W-shape pressure drop derivative curve is flat and small. The two simulated derivative curves of pressure drop for the cases without secondary fractures and without capillary imbibition deviate upward than that of the basic model, indicating that the oil-water displacement effect becomes worse in these two cases under the same pressure drop condition.

The established shut-in pressure drop model was used to fit the history of five typical fractured horizontal wells (J20, J21, J22, J23 and J24) in Jimsar shale oil reservoir of Xinjiang. The average fracture parameters of main hydraulic fractures, secondary fractures and reservoir, and the volumes of oil-water displacement were obtained from inversion by adjusting parameters.

Wells J20-J24 were all treated by hydraulic fracturing, with fracturing-fluid used of 34 408-63 589 m³ and fracturing cycle of 11-14 d. The wells were shut in for 29-57 d after fracturing. The wellhead tubing pressure of each well was continuously monitored during the shut-in period, as shown in Fig. 7a, Fig. 8a, Fig. 9a, Fig. 10a and Fig. 11a. A rapid pressure drop was observed at the early stage of shut-in. The pressure drop was mainly in the first day of shut-in, after which the pressure drop slowed down, which is consistent with the simulation results of this model. The fitted pressure drop curves of the five wells are shown in Fig. 7b, Fig. 8b, Fig. 9b, Fig. 10b and Fig. 11b. The pressure-drop derivative curves of the five wells are in a W-shape, which is consistent with the simulation results. According to the parameters from the fitting of pressure drop (Table 1), the main hydraulic fractures were 100-126 m in half-length and 5-10 μm²·cm in conductivity, the densities of secondary fractures were 0.42-0.93 fractures/m^-2, the matrix permeability was (0.005-0.012)×10^-3 μm², the volume of water imbibed by matrix was 4859-15 648 m³, and volume of oil replaced into fractures was 687-8013 m³.

The interpreted results of wells J20 and J23 were compared with micro-seismic monitoring records and numerical simulation results of geological engineering coupled with hydraulic fracturing treatment. According to micro-seismic monitoring, the hydraulic fractures of J20 and J23 horizontal wells are 116.5 m and 145 m in average half-length, respectively, 130.75 m on average. From the simulation of hydraulic fracturing, the hydraulic fractures of the two wells have an average half-length of 123.5 m^[30]. The half-length of main hydraulic fractures: 100-126 m, interpreted by this model are close to the micro-seismic monitoring records and numerical simulation results of geological engineering coupled with hydraulic fracturing treatment. Because fracture closure occurs during shut-in, the values from this model are slightly smaller than those from the other methods, which indicates that the fracture parameters obtained from this model shows the authentic effect of hydraulic fracturing.

According to the proposed shut-in pressure drop model coupling wellbore-fracture network-matrix system and numerical simulation results, the conclusions from this study are drawn as follows. In the initial stage of shut-in, after-flow dominates, with the maximum flow rate. In the early stage of shut-in, the channeling and leak-off in the fracture system dominate, and the fractures tend to close gradually. The middle stage of shut-in is dominated by matrix imbibition and oil replacement, when the leak-off of fracturing fluid and oil replacement from matrix tend to reach equilibrium. In the late stage of shut-in, the pressure wave reaches the reservoir boundary, and the leak-off rate of fracturing-fluid and oil replacement rate decrease simultaneously till to zero.

By analyzing the displacement rates of water and oil between different media during shut-in, the pressure drop characteristic curve is divided into eight stages: the stage dominated by wellbore after-flow, the stage dominated by channeling in fractures, the stage dominated by the leak-off from the fracture system, the stage of the fracture system storage, the stage of matrix imbibition and oil replacement, the stage of the steady increase of water-oil replacement efficiency, the stage of the decline of water-oil replacement efficiency, and the stage of water-oil replacement equilibrium.

Secondary fractures and capillary imbibition affect the shape of pressure drop derivative curve. The W shape does not appear without secondary fractures, while the pressure drop derivative curve for oil-wet formation without capillary imbibition shows nonstandard W-shape, with the second valley flat and small. Secondary fractures and imbibition related to wettability are the dominant effects of water-oil replacement.

The established shut-in pressure drop model was used to fit the pressure drop records from the oil field to obtain the parameters of main hydraulic fractures, secondary fractures and reservoir, and the volumes of oil-water replacement. The comparison of the half-lengths of the main hydraulic fractures from this model with the micro-seismic monitoring records and numerical simulation results of geological engineering coupled with hydraulic fracturing treatment shows that the fracture parameters obtained from this model shows the authentic effect of hydraulic fracturing.

C_F, C_f, C_m—the porosity compressibility coefficients of main hydraulic fracture, secondary fracture and matrix, Pa^-1;

d_F, d_f, d_m—the stress sensitivity factors of permeability of main hydraulic fracture, secondary fracture and matrix, Pa^-1;

D—the longitudinal migration distance of fluid, m;

g—the gravitational acceleration, m/s²;

K_F, K_f, K_m—the absolute permeability of main hydraulic fracture, secondary fracture and matrix, m²;

K_0,F, K_0,f, K_0,m—the initial permeability of main hydraulic fracture, secondary fracture and matrix, m²;

K_ro,F, K_ro,f, K_ro,m—the oil-phase relative permeability of main hydraulic fracture, secondary fracture and matrix, non-dimensional;

K_rw,F, K_rw,f, K_rw,m—the water-phase relative permeability of main hydraulic fracture, secondary fracture and matrix, non-dimensional;

p₀—original reservoir pressure, Pa;

p_c,m—the capillary pressure of matrix, Pa;

p_o,F, p_o,f, p_o,m—the pressures of oil-phase in main hydraulic fracture, secondary fracture and matrix, Pa;

p_w,F, p_w,f, p_w,m—the pressures of water-phase in main hydraulic fracture, secondary fracture and matrix, Pa;

p_wi,F, p_wi,f, p_wi,m—the pressures of water-phase in main hydraulic fracture, secondary fracture and matrix at the beginning of well shut-in, Pa;

p_wf—bottom-hole flow pressure, Pa;

q_o—the oil flux between wellbore and main hydraulic fracture, kg/(m³·s);

q_w—the injection rate of water-phase, kg/(m³·s);

q_w,WF, q_w,Ff, q_w,Fm, q_w,fm—the water fluxes between wellbore and main hydraulic fracture, between main hydraulic fracture and secondary fracture, between main hydraulic fracture and matrix, and between secondary fracture and matrix, kg/(m³·s);

S_o,F, S_o,f, S_o,m—the oil saturations of main hydraulic fracture, secondary fracture and matrix, %;

S_w,F, S_w,f, S_w,m—the water saturations of main hydraulic fracture, secondary fracture and matrix, %;

S_wi,F, S_wi,f, S_wi,m—the water saturations of main hydraulic fracture, secondary fracture and matrix at the beginning of well shut-in, %;

t—time, s;

x, y, z—the rectangular coordinate system, m;

α₁, α₂, α₃, α₄—the shape factors between wellbore and main hydraulic fracture, between main hydraulic fracture and secondary fracture, between main hydraulic fracture and matrix, and between secondary fracture and matrix, m^-2;

η_o, η_w—the viscosities of oil-phase and water-phase, Pa·s;

ρ_o, ρ_w—the densities of oil-phase and water-phase, kg/m³;

ϕ_F, ϕ_f, ϕ_m—the porosities of main hydraulic fracture, secondary fracture and matrix, %;

ϕ_0,F, ϕ_0,f, ϕ_0,m—the original porosities of main hydraulic fracture, secondary fracture and matrix, %.