The water–oil exchange phenomenon in hydraulically fractured shale reservoirs involves two-phase fluids, that is, oil and water, and multiple interconnected media, that is, main hydraulic fractures, induced natural fractures, and matrix pores. In the treatment of hydraulic fracturing, a large amount of water is injected into the wellbore, through the perforating hole into the reservoir to open fractures. When the well is shut-in after fracturing, the wellhead will not be injected anymore and the bottom-hole pressure will diffuse in the reservoir. At the same time, 100% water-saturated hydraulic fractures with high-pressure discharge water into the matrix through induced natural fractures and almost oil-statured shale matrix pores suck water, that is, the water imbibition under the potential difference, which includes hydraulic, capillary, and osmotic pressures (Wang et al., 2017b). Simultaneously, the oil in matrix pores is replaced to hydraulic fractures as compensation. Figure 1 shows a sketch of the oil–water exchange process.

A multi-pore medium comprising of main hydraulic fractures (F), induced natural fractures (f), and matrix pores (m) is proposed for characterizing hydraulically fractured shale reservoirs. The grid representation is exhibited in Figure 2. As shown in the figure, the whole medium is separated into two layers, with the lower one representing m and the upper one representing f. A shape factor (Kazemi et al., 1992) is used to quantify the complexity of f. The refined grids with high conductivity in the upper layer are used to characterize F, which is an identical transverse fracture with height, width, and length. F directly connects to the wellbore (W). In the soaking process, water suck in and oil discharge occurs between two adjacent porous media, as shown in Figure 2. q _w ^Ff and q _w ^fm represent the water exchange between F and f and between f and m, respectively, while q _o ^fF and q _o ^mf represent the oil exchange between F and f and between f and m, respectively. The water–oil exchange in this triple-porosity system is a hydrodynamic equilibrium process. It takes time for oil breakthrough in hydraulic fractures (t ^F _o,brh ) and water sucking termination in the matrix (t ^m _w,eqm ).

Following equations describe fluid mass transfer in F, f, and m, respectively. The subscript j represents water (w) and oil (o). F ∂ ( ρ j ϕ F S j F ) ∂ t = − ∇ ( ρ j v j F ) + q j W F − q j F − f , (1) where ρ j is density [kg/m³]; ϕ F is the porosity of F [-]; S j F is the fluid saturation in F [-]; v j F is the velocity in F [m/s]; q j W F is the fluid source sink term [kg/m²/s]; and q j F − f is the fluid transfer term linking F and f, and the flow direction is different between w and o [kg/m²/s]. v j F = − k F k r j F η j ∇ p j F , (2) where k F is absolute permeability of F [m²]; k r j F is fluid relative permeability of F [-]; p j F is hydraulic pressure in F [Pa]; and η j is viscosity [Pa·s]. q j F − f = α 1 ρ w k F k r j F h η j ( p j F − p j f ) , (3) where α 1 is the shape factor between F and f, and the expression and derivation of α 1 is presented in Supplementary Appendix SA [m⁻²] and p j f is hydraulic pressure in f, and when j = w, capillary pressure is considered [Pa]. q j W F = α 3 ρ w k F k r j F h η j B j ( p w f − p j F ) , (4) where α 3 is the shape factor which links F and W, and the expression and derivation of α 3 is presented in Supplementary Appendix SA [m⁻²] and B j is the fluid formation volume factor [-]; and p w f is the following pressure of bottom-hole [Pa]. h ∂ ( ρ j ϕ f S j f ) ∂ t = − ∇ ( h ρ j v j f ) + q j F − f − q j f − m , (5) where ϕ f is the porosity of f [-]; S j f is the fluid saturation in f [-]; v j f is the velocity in f [m/s]; and q j f − m is the fluid transfer term linking f and m, and the flow direction is different between w and o [kg/m²/s]. v j f = − k f k r j f η j ∇ p j f , (6) where k f is the absolute permeability of f [m²] and k r j f is the fluid relative permeability of f [-]. q j f − m = α 2 ρ j k f k r j f h η j ( p j f − p j m + p π ) , (7) where α 2 is the shape factor linking f and m, and the expression of α 2 is presented in Supplementary Appendix SA [m⁻²]; p j m is the hydraulic pressure in m, and when j = w, capillary pressure is considered [Pa]; and p π is the osmotic pressure only for w [Pa]. p π = λ R T V W ln x f x m , (8) where V w is the molar volume of water [m³/mol]; T is temperature [K]; R is ideal gas constant [J/(mol·K)]; λ is membrane efficiency [-]; x f is the water activity of f [-]; and x m is the water activity of m [-] m h ∂ ( ρ j ϕ m S j m ) ∂ t = − ∇ ( h ρ j v j m ) + q j f − m , (9) where ϕ m is the porosity of m [-]; S j m is the fluid saturation in m [-]; and v j m is the velocity in m [m/s]. v j m = − k m k r j m η j ∇ ( p j m − p π m ) , (10) where k m is the absolute permeability of m [m²] and k r j m is the fluid relative permeability of m [-].

Closed boundary is designed for outer boundary, while initial condition is set according to the end of injection process. The mathematical model is solved with the use of the finite difference method with the detailed algorithm presented in Supplementary Appendix SB. After the model solution, four key variables are obtained, that is, q _w ^Ff , q _w ^fm , q _o ^mf , and q _o ^fF , after which oil–water exchange efficiency and volume dynamics with soaking time can be calculated.

The oil–water exchange efficiency, which is the water imbibition flux divided by the oil drainage flux, can be calculated below.

Oil drainage efficiency of the fracture is as follows: E t s F = q o f F / q w F f , (11) where q _o ^fF is the oil exchange between F and f and q _w ^Ff is the water exchange between F and f.

Oil drainage efficiency of the matrix is as follows: E t s m = q o m f / q w f m , (12) where q _o ^mf is the oil exchange between f and m and q _w ^fm is the water exchange between f and m.

The oil–water exchange volume, which is the cumulative flux rate in the soaking period, can be calculated below:

Water sucking volume into the matrix is as follows: Q w m = ∑ q w f m ⋅ t s , (13) where t _s is time.

Oil drainage volume into hydraulic fractures is as follows: Q o F = ∑ q o f F ⋅ t s , (14)

Oil drainage volume into natural fractures is as follows: Q o f = ∑ q o m f ⋅ t s − ∑ q o f F ⋅ t s . (15)

Three key time points are also determined. The oil–water exchange equilibrium time (t ^m _w,eqm ) is the time when q w f m = 0 . The steady exchange efficiency time (t ^F _ow,std ) is the time when E t s F = 1 . The time of oil breakthrough in fracture (t ^F _o,brh ) is the time when q o f F > 0 .

Numerical simulation of the oil–water exchange is conducted with the use of this model. The length of the horizontal well is 1200 m, the total fracture stages are 15, and the well lies in the reservoir center (1500 m × 560 m×40 m). In every stage, there are four identical transverse hydraulic fractures with a fracture half-length of 140 m and a fracture spacing of 20 m along the horizontal wellbore.

The relative permeability and capillary force (Brooks and Corey, 1964) are used to represent the formation wettability, so the curves of each domain are set on the basis of the typical water-wet formation in the Bakken shale oil basin (Fakcharoenphol et al., 2014). The porosity and permeability of the fracture system are assigned according to Aguilera chart of compressibility coefficients (Aguilera, 1999). Other parameters input are detailed in Table 1. The fracturing fluid pumping scheme is simulated as the injection of water with pressure-dependent fracture porosity curve increases. The results of saturation as well as fluid pressure are input as initial conditions for the following 1 year soaking simulation.

Figures 3, 4 display the evolution of pressure profiles during the soaking periods in the fracture and the matrix, respectively. Results shows that during the well soaking, the water phase pressure inside and near the hydraulic fracture declines, while the distant pressure in the natural fracture and the matrix increases. That is a typical bottom-hole pressure diffusion and reservoir energy storage phenomena in soaking periods.

Figures 5, 6 display the evolution of saturation profiles during the soaking periods in the fracture and the matrix, respectively. It shows that during the well-soaking process, the water saturation in the fracture decreases, while the distant water saturation in the matrix increases continually. That is a typical water imbibition phenomenon in soaking periods.

Figures 7, 8 display the water and oil fluxes among F, f, and m, respectively, during soaking periods. It indicates that F discharges water to f, while m sucks water from f and discharges oil into the fracture system. The end time of F discharging water is about the 110^th day and the imbibition equilibrium time of m is about 169 days. The oil breakthrough time is about 1.29 days.

As shown in Figure 9, the simulated oil drainage efficiency of the fracture system shows a hump trend in the first several weeks, while the exchange efficiency of the matrix system shows a slightly downward trend. The maximum values of exchange efficiency occur at the end time of F discharging water and m sucking water, respectively. The steady exchange efficiency of F occurs at the soaking time of 17.25 days and lasts to 64 days, during which the value of exchange efficiency is 0.805. Figure 10 shows the evolution of water–oil exchange volumes with soaking time. It indicates that the three volumes of water sucking into m, oil drainage into F, and f all increase with soaking time, but the increment tends to slow down from the 64^th day on. At the end of 1 year soaking, there are 1486 m³ of water sucking into the matrix and 1056 m³ of oil being replaced into the fracture system (21 m³ in F and 1035 m³ in f).

Figure 11A displays the simulated bottom-hole flowing pressure dynamics of the base case during 1 year soaking. To describe the water–oil flow dynamics, pressure falloff analysis theory is applied. Figure 11B displays the calculated pressure falloff and Bourdet derivatives. The Bourdet derivatives show a W-shaped curve. The first derivative rising stage in W-shaped derivative curve starts from 1.29 days, corresponding to the simulated oil breakthrough time t ^F _o,brh. In the later stage, the fluctuations of the derivative curve, like the inflection points at 5.6 and 8.7 days, are consistent with the trend change points on the oil drainage efficiency curve of the hydraulic fracture, that is, the blue curve shown in Figure 9. The steady exchange efficiency of F, which occurs from 17.25 to 64 days, is reflected by a constant slope of the second derivative rising stage in the W-shaped derivative curve. After soaking for 64 days, the pressure drop is very small, which exhibits the derivative curve falling. Until soaking 110 days, there is no pressure drop, and the derivative curve drops to zero.

The fluid mass transfer between different media and the pressure of these media are difficult to monitor, because monitoring techniques in the soaking period are limited. Therefore, the wellhead or bottom-hole pressure falloff is the only index that can be used to validate the model.

Shale formation C is characterized by the following reservoir properties: initial pressure = 38 MPa; effective thickness = 10∼40 m; permeability = 0.001∼0.01 md; porosity = 0.07∼0.11; and induced natural fracture density = 0.16∼0.25 m⁻². The induced natural fracture density is obtained from the field report. It is converted to the α ₂ in the model (α ₂ = 2.5∼5 m⁻²) based on previous studies (Kazemi et al., 1992; Yan and Mi, 2017). This formation is stimulated by hydraulic fracturing treatment for commercial oil production. The volume of fracturing fluid injected to each well ranges from 7,920 to 6,3589 m³; the volume of proppant added to each well ranges from 440 to 4,550 m³; the fracture stages for each well range from 5 to 36. The wells in this formation are shut-in from 19 to 61 days after stimulation. During the well soaking period, the wellhead pressure is monitored continuously. The wellhead pressure records of 10 wells during post-fracturing soaking periods are shown in Figure 12A. The Bourdet derivatives of recorded wellhead pressure falloff, shown in Figure 12B, exhibit the typical W-shaped curves, which are consistent with our simulated results in Figure 11B. Based on the simulated relevance of time between the oil–water exchange and the pressure falloff derivatives above, it can be inferred that oil breakthrough already happens for most of wells since the first derivative rising stage in their W-shaped derivative curves exists from the initial 3 days recorded pressure data. However, the derivative falling after the second rising stage in their W-shaped derivative curves is not prevalent, indicating that water–oil exchange efficiency of well is still in a rising stage and it is still early to open the well for production.