Numerical Investigation of Oil–Water Exchange Behaviors in Shale During Post-Fracturing Soaking Periods

Fracturing fluid imbibition and retention are treated as a main mechanism for oil production from shale reservoirs. However, the oil–water exchange phenomenon during post-fracturing soaking periods has not been thoroughly studied. In this study, a water–oil flow model is built to investigate the water imbibition and oil drainage phenomenon in hydraulically fractured shale. With the developed numerical simulator, the main characteristics of post-fracturing soaking, that is, pressure diffusion, water imbibition, and especially, the oil–water exchange behavior are simulated. Three key time points, that is, oil–water exchange equilibrium, steady exchange efficiency, and oil breakthrough in fracture are found. The oil–water exchange efficiency and exchange volume are also calculated. Moreover, the proposed model is validated by field wellhead pressure dynamics, indicating a relevance of time between the oil–water exchange efficiency and the wellhead pressure falloff derivatives. Finally, the influences of shale permeability, wettability, fracture complexity, and oil viscosity on the oil–water exchange behavior are investigated. Results indicate that the matrix of oil-wet shale almost does not suck water and discharge oil, and only the oil in natural fractures exchanges with the water in hydraulic fractures. The water-wet shale with low permeability, high oil viscosity, and few natural fractures needs extra soaking time to achieve good oil–water exchange performance. The suitable soaking period for the water-wet base case in this study is from 17.25 to 169 days, among which 64 days is the optimal soaking time.


INTRODUCTION
Compared with conventional oil reservoirs, shale has relatively high clay content and covers a wide range from 16.8 to 70.1% (Yang et al., 2013). The development of shale oil reservoirs heavily relies on multistage hydraulic fracturing technology (Zhou et al., 2019). After the hydraulic fracturing, the initial production of fractured shale oil wells is good, but the stable production period is quite short (Zou et al., 2020). For pursuing a high oil production rate, a long-term well shut-in has gradually been adopted as an effective mode for field practice, instead of flowing back immediately after hydraulic fracturing because it is believed that fracturing fluid imbibition may strengthen and oil-water exchange may happen during the soaking period. Although extended well soaking may intensify fracturing fluid retention, which is proved by field practice and laboratory experiments. After all, sacrificing water recovery for more oil is acceptable.
The investigation of imbibition behavior is mainly based on laboratory experiments and a pore/core scale analysis. In general, fracturing fluid imbibition behavior can be related to many factors, including rock type (Xiong et al., 2013a;Xiong et al., 2013b;Ren et al., 2015), mineral composition (Akin and Kovscek, 1999;Ali Habibi et al., 2015;Ge et al., 2015), pore structure (Liu and Dai, 2008), formation wettability (ZhuJu et al., 2002;Wang et al., 2012), formation temperature (Li et al., 2011;Zhou et al., 2014), the bedding structure (Ren et al., 2015), connate water saturation (Gao and Hu, 2016), formation fluid (Wang, 2000;Wang and Zhang, 2000), and fracturing fluid type (Paktinat et al., 2005;Penny et al., 2005;Li, 2006;Roychaudhuri et al., 2013). Based on the results of experiments, formation wettability is the most important factor for fluid imbibition since it determines the direction of imbibition, and the extent of wettability determines the oil recovery. The contact area affects the imbibition speed, that is, the larger the contact area is, the faster the imbibition speed is (ZhuJu et al., 2002). Moreover, fluid properties are critical to the imbibition behavior. The content of polar substance in oil may alter the formation wettability and further change the imbibition behavior (Wang, 2000;Wang and Zhang, 2000). In practice, surfactants are added into the fracturing fluid. The existence of surfactant can reduce the interfacial tension (Li, 2006;Roychaudhuri et al., 2013), making the rock more hydrophilic to enhance the capacity of water suction, while microemulsion as a cleanup additive can make the rock wettability to water from strong to weak and then inhibit the capacity of water suction (Paktinat et al., 2005;Penny et al., 2005). Besides, the porosity and permeability of formation are very important to fluid imbibition. It shows that the capillary pressure increases with the decrease of porosity and permeability, which makes the driving force of imbibition to increase (Liu and Dai, 2008).
Besides the experiment-based pore/core scale analysis, the well/reservoir scale imbibition, and the well shut-in simulation mainly focuses on shale gas reservoirs. Wang et al. (2017a) studied the phenomenon of fluid imbibition under the driving forces of osmosis and capillarity and investigated the influence of imbibition on the microfractures induced by hydraulic fracturing in shale gas reservoirs. Their research result shows that clay minerals have the function of a semipermeable membrane and the capacity of water suction is stronger than organic matter and other minerals. Fakcharoenphal et al. (2013) shows that fracturing fluid imbibition in shale can generate microfractures, and it is related with the physical characteristics of reservoir itself and the composition of fracturing fluids. Roychaudhuri et al. (2013) shows that shale reservoir has the characteristics of mixed wetting, which has an important impact on the fracturing fluid imbibition. Odumabo et al. (2014) shows the relation between imbibition volume and soaking time. It shows that the change of permeability in the invasion area of the fracture surface and the distribution of water saturation near the fracture surface are two important factors affecting the post-fracturing shale gas production. Fakcharoenphol et al. (2014) studied the effect of salinity on the fluid imbibition. It shows that osmotic pressure is an important driving force for fluid imbibition. Wang and Leung (2015) established a triple-porosity two-phase flowback model to study the control mechanism of fracturing fluid retention. Their simulation results show that the initial production rate can be increased by extension of soaking time, but for a long-term production rate, soaking is not helpful. Zhang et al. (2017) established a triple-porosity numerical model for hydraulically fractured shale gas wells, which considers the effects of imbibition, stress sensitivity, and gravity differentiation. Their simulation results show that the spontaneous imbibition rate of fracturing fluids in a tight matrix is very low during the well-soaking process. Although shale reservoir has high capillary force, the fracturing fluid still mainly remains in the fracture and near the fracture area after 100 days of soaking. Their simulation results show that although the long-term soaking increases the initial gas production, the long-term cumulative gas production decreases with the extension of soaking time.
From the perspective of modeling and simulation, water imbibition and oil drainage has not been systematically coupled with different driving mechanisms, such as hydraulic pressure difference, wettability-dominated capillarity, and chemical osmosis. To investigate the oil-water exchange behavior in a well/reservoir scale, we build a fracture-matrix water-oil flow model to simulate the fracturing fluid imbibition and the oil drainage process during the post-fracturing soaking periods. We try to use the numerical simulation results to investigate the efficiency and equilibrium time of oil-water exchange under different shale permeability, wettability, fracture complexity, and oil viscosity conditions, which help to optimize the well-soaking time for hydraulically fractured shale oil wells.

Mechanism Description
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 highpressure 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.

Mechanism Model Development
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 ).

MATHEMATICAL SOLUTION
Assumptions made for this study are as follows. 1) Sodium chloride is the only solute and totally dissolves in either formation brine or fracturing fluid; 2) gravity is not considered; 3) mass transfer occurs inside F, f, m, and between the adjacent media; and 4) the system is isothermal. According to the aforementioned physical model and assumptions, a mathematical model is developed as follows.

Oil-Water Flow Model
Following equations describe fluid mass transfer in F, f, and m, respectively. The subscript j represents water (w) and oil (o). F where   linking F and f, and the flow direction is different between w and o where α 1 is the shape factor between F and f, and the expression and derivation of α 1 is presented in Supplementary Appendix SA [m −2 ] and p f j is hydraulic pressure in f, and when j w, capillary pressure is considered [Pa].
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 −2 ] andB j is the fluid formation volume factor [-]; and p wf is the following pressure of bottom-hole [Pa].
where k f is the absolute permeability of f [m 2 ] and k f rj is the fluid relative permeability of f [-].
whereα 2 is the shape factor linking f and m, and the expression of α 2 is presented in Supplementary Appendix SA [m −2 ]; p m j 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].
where k m is the absolute permeability of m [m 2 ] and k m rj is the fluid relative permeability of m [-].

Model Coupling and Variables Solution
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: 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: 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: where t s is time.
Oil drainage volume into hydraulic fractures is as follows: Oil drainage volume into natural fractures is as follows: Three key time points are also determined. The oil- 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.

Simulation Results of Pressure and Saturation Fields
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 bottomhole 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.  (Fritz and Marine, 1983;Cheng et al., 2009;Fakcharoenphol et al., 2014;Wang et al., 2020   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 3 of water sucking into the matrix and 1056 m 3 of oil being replaced into the fracture system (21 m 3 in F and 1035 m 3 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.

Model Validation
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 −2 . The induced natural fracture density is obtained from the field report. It is converted to the α 2 in the model (α 2 2.5∼5 m −2 ) 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 3 ; the volume of proppant added to each well ranges from 440 to 4,550 m 3 ; 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.

SENSITIVITY ANALYSIS
The influences of shale permeability (k m ), wettability, fracture complexity (α 2 ), and oil viscosity (µ o ) on the oil-water exchange behavior are investigated. The evaluation indicators include oil breakthrough time (t F o,brh ), steady exchange efficiency time (t F ow,std ), exchange equilibrium time (t m w,eqm ), and water-oil exchange efficiency volumes (Q F o, Q f o , and Q m w ). The oil-water relative permeability and capillary pressure curves for the sensitivity analysis of oil-wet shale formation are assigned below (Wang, 2000;Li, 2006). Figure 13 shows the times of imbibition equilibrium, oil breakthrough, and steady exchange efficiency with different shale prosperities during well soaking. It indicates that µ o increasing from 0.36 to 1.62 causes t F o,brh , t F ow,std , and t m w,eqm to increase, while both α 2 increasing from 0.03 to 10 and k m increasing from 0.001 to 0.01 cause t F o,brh , t F ow,std , and t m w,eqm to decrease. By contrast to the water-wet base case, the oil-wet case shows a shorter imbibition equilibrium time (t m w,eqm 0.25 days) and a longer oil breakthrough time (t F o,brh 1.75 days). The steady exchange efficiency for the oil-wet case is about 0.01, which is far below that of the water-wet base case. Figure 14 shows the 1-year oil-water exchange volumes with different shale prosperities. The simulation results indicate that µ o  to increase first and then decrease, which indicates a complicated imbibition-replacement behavior in the fracture-matrix system. By contrast to the water-wet base case, the oil-wet case shows that the matrix sucks water 36.72 m 3 within the first 0.25 days and then discharges water for 1 year, which leads to a negative watersucking volume (Q w m −851 m 3 ) and a negative oil discharge volume (Q o f −697 m 3 ). The oil drainage volume in the hydraulic fracture is very small: Q o F 0.003 m 3 . Based on our simulation results above, the oil-water exchange behavior during soaking periods can be summarized and used to direct the engineering practice for determination of optimal wellsoaking time. For the water-wet base case, the water imbibition and oil drainage occur in the matrix at the beginning of the soaking periods, which exhibits that the water in F enters m through f, but the oil breakthrough in F will be delayed to 1.29 days. The water sucking equilibrium of the matrix occurs on the 169 th day, which means that the oil-water displacement driven by hydraulic pressure difference ends, after which the water flow into the matrix is driven by spontaneous imbibition. The imbibition equilibrium time of the matrix, that is, 169 days, can be the upper limit soaking time. Before that, the water-oil exchange efficiency shows a hump trend. The imbibition efficiency tends to be stable from 17.25 to 64 days, and the stable value is 0.805. The volume of water-oil exchange from 64 to 169 days is very small, which is mainly the amount that the oil already discharged from f and enters F. So the suitable soaking period for the water-wet base case is from 17.25 to 169 days, among them 64 days are the optimal soaking time.
From the sensitivity simulation, it can be inferred that the water-wet shale reservoir with low permeability, high oil viscosity, and few natural fractures needs extra soaking time to achieve good oil-water exchange performance. However, for the oil-wet shale reservoir, the matrix almost does not suck water and   Frontiers in Earth Science | www.frontiersin.org October 2021 | Volume 9 | Article 735972 8 discharge oil, but only the oil in natural fractures exchanges with the water in hydraulic fractures. So the water-oil exchange efficiency is very low, and extending the soaking time cannot achieve good oil-water exchange performance.

CONCLUSION
In this study, a water-oil flow mechanism model is built to investigate the water-oil exchange behaviors in hydraulically fractured shale during soaking periods. Some understandings of the key points are summarized below.
1) The simulation results prove the occurrence conditions for water-oil exchange, that is, water-wet rock can make the oil in the matrix replaced by the water in fractures, while for oil-wet rock, the oil in the matrix cannot be replaced out by well soaking and the oil in natural fractures will be further squeezed into matrix pores with the extension of soaking time, resulting in production difficulty. 2) The simulation results reflect the water-oil exchange behavior in water-wet shale, that is, the water imbibition and oil drainage occur in the matrix at the beginning of the soaking periods, but the oil breakthrough will be delayed to 1.29 days. The water-sucking equilibrium in the matrix occurs on the 169 th day, which means that the oil-water exchange driven by hydraulic pressure difference ends, after which the water flow into the matrix is driven by spontaneous imbibition.
3) The simulation results indicate that low matrix permeability and high oil viscosity are detrimental to oil-water exchange, which exhibits low exchange efficiency, late oil breakthrough, and longer equilibrium time. On the contrary, the existence of natural fractures is advantageous to oil-water exchange, which is characterized by short imbibition equilibrium time, fast oil breakthrough, and high oil exchange efficiency. 4) Simulation results guide the determination of optimal soaking time. The steady exchange efficiency time (t F ow,std ) and imbibition equilibrium time (t m w,eqm ) can be the lower and upper limit soaking time, respectively. Among this period, the optimal soaking time can be further determined according to the expected oil discharge rate of the matrix and oil drainage rate of hydraulic fractures. 5) The model simulation results also indicate a relevance of time between the oil-water exchange efficiency of hydraulic fractures and the pressure falloff derivatives. This indication can be used for determining oil-water exchange efficiency of actual wells.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material; further inquiries can be directed to the corresponding author.