Method for Calculating Productivity of Water Imbibition Based on Volume Fracturing Stimulations of Low Permeability Reservoirs

Spontaneous water imbibition is an important mechanism in water-wet fractured reservoirs. For volume-fractured reservoirs, to evaluate the oil productivity and oil recovery through water counter-current imbibition, we propose an analytical method for optimizing the reservoir volume fracturing scheme. Based on the two-phase ﬂ uid ﬂ ow di ﬀ erential equation for capillary force, a three-dimensional water imbibition productivity equation is derived analytically. The equation for the water imbibition productivity considering the fracture network is obtained. A numerical model is constructed to verify the validity of the average capillary di ﬀ usivity coe ﬃ cient and the results of the analytical model. By applying this method to a low permeability reservoir, after volume fracturing and water ﬂ ooding hu ﬀ and pu ﬀ , the relationship between the tenth year ’ s oil recovery and oil production rate and the length, width, and density of the fracture network is predicted, which gives an optimization of the ﬁ eld fracturing construction scale. The results show that the length and width of the fracture network should be no less than 50% of the well spacing and row spacing to obtain a reasonable production. Considering the fracturing technique and economic feasibility, the higher the density of the fracture network, the better the production obtained. Through hydraulic volume fracturing and water ﬂ ooding hu ﬀ and pu ﬀ , water imbibition is brought into full play and the 10 year oil recovery is increased by 6% – 8% in this area.


Introduction
The volume fracturing technique is a very useful tool for developing low permeability and unconventional oil and gas reservoirs. It can significantly increase the productivity of wells and improve the final oil recovery. The volume fracturing stimulation mechanism involves forming a fracture network interwoven with the main fractures and multistage secondary fractures, which maximizes the stimulated reservoir volume and effectively reduces the fluid seepage distance from the pores to the fractures [1][2][3][4][5][6]. In water-wet fractured reservoirs, the difference in the capillary pressures of the matrix media and the fractures provides the main driving force for water imbibition (Figure 1), and this has become an efficient oil recovery mechanism [7][8][9][10][11][12].
According to the flow direction of the oil and water, imbibition can be divided into two modes: countercurrent and cocurrent flows. Intuitively, countercurrent imbibition is defined as the case when water is imbibed into the rock while oil is expelled in the opposite direction of the water flow. In contrast, cocurrent imbibition is defined as the case when oil flows downstream of the water front and in the same direction as the water flow. A large number of studies [13][14][15][16][17] have shown that when the block is completely immersed in water or is surrounded by water in fractures, countercurrent flow is the dominant mechanism of the spontaneous water imbibition.
Over the years, several analytical models have been proposed to quantify the oil recovery in the water imbibition progress [18][19][20][21][22][23][24][25][26][27]. Aronofsky et al. [18] proposed an empirical model, in which the oil recovery is exponentially related to the imbibition time. Mattax and Kyte [19] defined a scaling group where the dimensionless time was used to scale the spontaneous imbibition. Later, Cuiec et al. [20][21][22][23] modified this model by introducing the viscosity ratio, characteristic length, and other parameters. Handy [24] reported that imbibition can be considered as a diffusion-like or piston-like process, and the volume of the imbibed water is proportional to the square root of the imbibition time. Zimmerman and Bodvarsson [25] and Li [26] used approximate analytical approaches to derive the fluid saturation distribution inside the matrix block. The Handy model has also been used in the semiquantitative analysis of shale gas [27]. Yang et al. [28] applied the imbibition rule and verified that the change in the salt concentration due to ion diffusion is proportional to the square root of time.
To analytically quantify the productivity of water imbibition, in this study, we begin by simplifying the imbibition diffusion equation to obtain the distribution of the water saturation in a three-dimensional core model. Then, the equation for calculating the productivity and oil recovery of water countercurrent imbibition is proposed. By generalizing the equation for water imbibition progress in volume fracturing reservoirs, we establish a quantitative relationship between the oil recovery factor and the network fracture parameters, such as the fracture density and fracture scales. Through this research, the oil recovery obtained using the water imbibition and the volume fracturing technique is quantitatively calculated, and a new method for optimizing the field volume fracturing scale is proposed.

Mathematical Model of Water Countercurrent
Imbibition. McWhorter et al. [29] first proposed the onedimensional incompressible flow equation for two-phase fluids considering capillary pressure: where S w is the water saturation, f w ′ ðS w Þ is the derivative of the water cut, Ф is the porosity; qðtÞ is the sum of the water and oil flow rates (m 3 ); A is the area of the contact surface of the matrix (m 2 ); x is the distance in the x direction (m); and t is the imbibition time (s). d w ðS w Þ is called the capillary diffusivity coefficient (m 2 /s), which is expressed as where K is the absolute permeability (10 -3 μm 2 ); k ro and k rw are the relative permeabilities of oil and water, respectively; μ w and μ o are the viscosities of water and oil (mPa·s), respectively; and P c is the capillary pressure of the oil and water phase (MPa). When countercurrent imbibition occurs, the flow rates of the two phases have the same value but opposite directions, so qðtÞ is equal to 0 in this case. The flow equation for countercurrent imbibition can be written as As can be seen from Equation (2), the capillary diffusivity coefficient is a function of the water saturation, and the relationship curve for d w and S w is shown in Figure 2. d w ðS w Þ is a bell-shaped function, and d w ðS wi Þ = d w ð1 − S or Þ = 0, where S wi is the initial water saturation and S or is the residual oil saturation of the core. Between these two endpoints, there is a certain saturation where d w has the highest value, d max .
Since Equation (3) is a nonlinear partial differential equation, the equation can be simplified into an ordinary differential equation and can be solved analytically only when the capillary diffusivity coefficient is a constant. When the capillary diffusivity coefficient is a variable, only approximate solutions can be obtained through numerical calculations. In this study, we used both numerical and analytical methods to solve the countercurrent imbibition equation.

Numerical
Method of Oil Productivity. First, the numerical method was applied to calculate the equation. The finite difference method was used for the numerical solution. We used D to represent d w ðS w Þ, which is a variable, and Equation (3) is rewritten as where S = ðS wm − S w Þ/ðS wm − S wi Þ, which is the normalized water saturation, and S wm = 1 − S or is the highest water saturation of the core. By extending Equation (4) into three-dimensional space, we obtained where D x , D y , and D z are the partial differentials of D in the x, y, and z directions, respectively. i, j, and k refer to the grid number in the x, y, and z directions, respectively; and n refers to the grid number of the time step. Using the finite 2 Geofluids difference method, the first order central difference of the variables x, y, and z were determined: In addition, the first-order forward difference was determined for the variable t: The five-point difference grid is shown in Figure 3, and Equation (5) becomes where Then, a one-dimensional mathematical model of countercurrent imbibition was established. As Figure 4 shows, the core is saturated in oil and immersed in water. Both the oil and water phases flow in the x direction through this core. The length of the matrix core is L, with the top, bottom, and right sides closed and only the left face open to the water. Gravity is ignored. At the initial time, the oil and water are in contact at position x = 0, where x = L is a sealed face. The equations are listed as follows:

Geofluids
For the 2-n-1 grid, the difference equation is where For the first grid, For the nth grid, To ensure the convergence of the iteration, that is, the algorithm's stability, according to the Von Neumann condition of the parabolic equation, For a given Δx, Δt needs to be small enough to satisfy Equation (15) in order to ensure the stability of the numerical calculation. Figure 5 and Table 1 show the relative permeability curves, capillary pressure curve, core properties, and fluid parameters. By encoding and running a MATLAB program, the distribution of the water saturation with respect to position x and time t through the core were obtained. The simulation results are shown in Figure 6, in which the x axis represents the relative distance from the original oil-water interface, and the y axis represents the imbibition time. At a certain time, as the relative distance increases, the water saturation decreases gradually. At a certain position, the water saturation increases with increasing imbibition time.
For a certain position and time during the imbibition process, we can obtain the value of the water saturation. Based on the distribution of the water saturation at each point, the average saturation of the matrix can be calculated, and the movable oil recovery factor of the water imbibition at a certain time can be obtained: where L is the length of the core (m) and S w is the average saturation of the core. The results of the numerical calculation of the movable oil recovery factor during the water imbibition progress are shown in Figure 6.

Analytical Method of Oil Productivity.
Owing to the large amount of calculations, the numerical method has a slow calculation speed and high convergence requirements, and thus, the analytical method is considered. However, the original problem needs to be solved, and Equation (3) can be written as an ordinary differential equation and an analytic solution can be obtained only when the capillary diffusivity coefficient is a constant. Thus, an average capillary diffusivity coefficient d is needed to represent the diffusion process.
The determination and verification of this average value is discussed later. When the capillary diffusivity 2.3.1. One-Dimensional Model. A one-dimensional model was constructed (Figure 7) as follows. We assumed that the oil and water phases flow along the 1D direction (x axis) and the matrix core has a length of L. The model is the same as the one built using the numerical method, except that both the left face and right face are open to allow contact between the oil and water phases. The sectional area of the contact surface is A. The water saturation of the matrix at the initial moment is S wi . The oilwater contact surface is located at x = 0 and x = L, where the water saturation is S wm . Thus, according to the 4 Geofluids hypothesis, the following mathematical model can be obtained: By separating the variables, the distribution law of the water saturation through the core can be obtained as follows: where t D = dt/ϕL 2 , x D = x/L are the dimensionless time and the dimensionless position, respectively. By integrating Equation (19), the average water saturation of the core can be obtained as follows: The oil recovery of water imbibition at t is The cumulative oil production QðtÞ can be calculated using Equation (21), and the derivative with respect to time can be taken to obtain the oil production: 2.3.2. Verification. The determination of the average capillary diffusivity coefficient d is very important, and the value of d should represent the imbibition process. By comparing the results of the numerical method and the analytical method using a certain constant of d, we found the best way to obtain the average value. In this way, the nonlinear dependency of the diffusion coefficient was linearized.
Taking the one-dimensional model as an example, we used the parameters in Figure 5 and Table 1 to calculate the oil recovery of this model, and the result is shown as the   Figure 8. Then, we determined the integral mean value of the capillary diffusion coefficient: The result obtained by calculating the oil recovery using Equation (21) is shown as the solid line in Figure 8. The comparison shows that except for the slight deviations in the initial and final stages of the imbibition process, the analytical results are in good agreement with the numerical results when the diffusion coefficient is taken as the integral mean value of d w ðS w Þ. The numerical method verifies the determination of the average diffusion coefficient. Therefore, all of the calculations below, which were conducted using the analytical method, are based on the mean integral value of the diffusion coefficient. Figure 9, a cuboid model was constructed. Both the oil and water flow along all three directions (x, y, and z axes). All of the six faces of the core are in contact with the water. L x , L y , and L z are the lengths of the core in the x, y, and z directions, respectively. The initial conditions are consistent with the one-dimensional model. Based on the equations for the one-dimensional model, the water saturation distribution, oil recovery, and production of the three-dimensional model were deduced.

Three-Dimensional Model. As is shown in
The water saturation distribution of the countercurrent imbibition in the three-dimensional model is

Geofluids
The oil recovery of the water imbibition at t is The oil production of the water imbibition is 2.4. Imbibition Oil Productivity considering the Facture Network. To predict the imbibition capacity of the fracture network ( Figure 10) formed by the volumetric fracturing in a low permeability reservoir, the following study was conducted based on a fracture-pore dual media model. The fracture parameters, such as the length, width, and density of the fracture network, were optimized to obtain a better oil productivity through water imbibition. Figure 11 shows the top view of this model. We discuss an oil well as an example. The hypotheses are as follows. The thickness of this reservoir is h, the lengths of the fracture network along the x axis and y axis are a and b, respectively. The numbers of fractures in the x-axis and y-axis directions are f x and f y , respectively. Since the fracture network divides the area near the wellbore into several matrix blocks, according to the number of contact surfaces with fractures, the matrix blocks are divided into three categories. The different types of matrix blocks correspond to the different methods for calculating the oil productivity (Table 2).
Based on the 3D oil productivity equation, we can obtain the oil recovery and production considering the fracture network.
The oil recovery of the water imbibition at t considering the volumetric fracturing is The oil production of the water imbibition considering the volumetric fracturing is

Results and Discussions
The case study focuses on a low permeability reservoir located in the Bohai Bay Basin. The original oil-in-place (OOIP) is about 17 MMBBL, the average porosity is 13%, and the permeability is 3:9 × 10 −3 μm 2 . The oil viscosity under the reservoir conditions is 2.3 mPa·s and the water vis-cosity is 0.5 mPa·s. This reservoir uses a rectangular well pattern with a well spacing of 300 m and a row spacing of 200 m. Based on the experimental data we determined that the residual oil saturation of this reservoir is 0.36 and the irreducible water saturation is 0.35. In addition, we used the integral mean value method to calculate the average capillary diffusivity coefficient (2:24 × 10 −8 m 2 /s). The average initial

Matrix block type Number
Two faces in contact with fractures (1 in Figure 10) 4 Three faces in contact with fractures (2 in Figure 10) Four faces in contact with fractures (3 in Figure 10) 8 Geofluids production of the wells was 63 BPD. The cumulative oil production during the 10-year evaluation period was about 1 MMBBL. The oil recovery factor of the water flooding during this period was 5.85%. Now the development mode has changed to volume fracturing with waterflooding huff and puff. One producing cycle is 1-year long, which includes 2 months of water injection, 1 month of well shut-in, and 9 months of production. We used Equations (27) and (28) to calculate the oil production and recovery of the water imbibition in order to design an optimized volume fracturing development plan. Figures 12 and 13 show the relationships between the length, width, and density of the fracture network and the oil recovery of the water imbibition. The results show that the oil recovery increases as the fracture density, length, and width of the fracture network increase. When the fracture network is longer than 150 m (half of the well spacing) or wider than 100 m (half of the row spacing), the degree of increase gradually flattens out. However, as the fracture density increases, the oil recovery continues to increases significantly. Figures 14 and 15 show the relationships between the oil production and the fracture length and width in log coordinates. In the initial stage, the oil production was high; however, as the imbibition time increased, the production decreased rapidly. When the length of the fracture network is less than 150 m or the width is less than 100 m, the   9 Geofluids production of the imbibition at the end of one cycle decreased to less than 13 BPD. The proposed optimized volume fracturing development plan is as follows. The length and width of the fracture network should be no less than 150 m (half of the well spacing) and 100 m (half of the row spacing), respectively; otherwise, the oil production and recovery factors will be significantly lower. Considering the volume fracturing technique and economic feasibility, the higher the density of the fracture network, the higher the pro-duction and the better the productivity of the water imbibition.

Conclusions
An analytical method was proposed to optimize the volume fracturing scheme of the reservoir and to evaluate the oil productivity and recovery by water countercurrent imbibition in volume fracturing reservoirs. To save calculation time and    provide a quantitative method for calculating the water imbibition productivity, the determination of the average capillary diffusivity coefficient d is very important and the value of d should represent the imbibition process. In this study, we determined the integral mean value of the capillary diffusion coefficient, and the accuracy of this analytical method was verified through numerical calculations.
For volume fracturing of reservoirs, the oil recovery is positively correlated with the length, width, and density of the fracture network, which allows for the optimization for the field fracturing construction scale. To achieve more production, the length and width should be no less than 50% of the well spacing and row spacing, respectively. Considering the fracturing technique and economic feasible, the higher the density of the fracture network, the better the production. Using hydraulic volume fracturing and cyclic water injection, water imbibition is brought into full play, and the tenth year's oil recovery is increased by 6%-8%.
Through this study, the laboratory scale research results for water imbibition were applied at the field volume fracturing scale, and in addition to a numerical method, we propose an approximate analytic method for optimizing the scale of the fracture network in volume fracturing reservoirs, which provides theoretical support for field volume fracturing.

Data Availability
The raw/processed data required to reproduce these findings cannot be shared at this time as the data also forms part of an ongoing study.

Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this study.