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

Authors: Lisha Zhao; Xin Li; Shuhong Wu; Zhongbao Wu; Min Tong
Journal: Geofluids
 (2021)

---

## Abstract

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 fluid flow differential 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 diffusivity coefficient and the results of the analytical model. By applying this method to a low permeability reservoir, after volume fracturing and waterflooding huff and puff, 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 field 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 waterflooding huff and puff, water imbibition is brought into full play and the 10 year oil recovery is increased by 6%–8% in this area.


---

## Body

## 1. 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–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–12].Figure 1
Water imbibition progress between the matrix and fracture.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–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–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–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.


## 2. Materials and Methods

### 2.1. Mathematical Model of Water Countercurrent Imbibition

McWhorter et al. [29] first proposed the one-dimensional incompressible flow equation for two-phase fluids considering capillary pressure:
(1)qtϕAfw′Sw∂Sw∂x−1ϕ∂∂xdwSw∂Sw∂x+∂Sw∂t=0,where Sw is the water saturation, fw′Sw is the derivative of the water cut, Ф is the porosity; qt is the sum of the water and oil flow rates (m3); A is the area of the contact surface of the matrix (m2); x is the distance in the x direction (m); and t is the imbibition time (s). dwSw is called the capillary diffusivity coefficient (m2/s), which is expressed as
(2)dwSw=−Kkrwkrokrwμo+kroμwdPcdSw,where K is the absolute permeability (10-3μm2); kro and krw are the relative permeabilities of oil and water, respectively; μw and μo are the viscosities of water and oil (mPa·s), respectively; and Pc 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 qt is equal to 0 in this case. The flow equation for countercurrent imbibition can be written as
(3)1ϕ∂∂xdwSw∂Sw∂x=∂Sw∂t.As can be seen from Equation (2), the capillary diffusivity coefficient is a function of the water saturation, and the relationship curve for dw and Sw is shown in Figure 2. dwSw is a bell-shaped function, and dwSwi=dw1−Sor=0, where Swi is the initial water saturation and Sor is the residual oil saturation of the core. Between these two endpoints, there is a certain saturation where dw has the highest value, dmax.Figure 2
Relationship curve of the capillary diffusivity coefficient and water saturation (the dotted line is the average capillary diffusivity coefficient’s value).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.


### 2.2. 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 usedD to represent dwSw, which is a variable, and Equation (3) is rewritten as
(4)∂∂xD∂S∂x=ϕ∂S∂t,where S=Swm−Sw/Swm−Swi, which is the normalized water saturation, and Swm=1−Sor is the highest water saturation of the core.By extending Equation (4) into three-dimensional space, we obtained
(5)∂∂xDx∂S∂x+∂∂yDy∂S∂y+∂∂zDz∂S∂z=ϕ∂S∂t,where Dx, Dy, and Dz 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 difference method, the first order central difference of the variables x, y, and z were determined:
(6)∂f∂x=fi+1−fi−12Δx,∂f∂y=fj+1−fj−12Δy,∂f∂z=fk+1−fk−12Δz.In addition, the first-order forward difference was determined for the variablet:
(7)∂f∂t=fn+1−fnΔt=δfΔt.The five-point difference grid is shown in Figure3, and Equation (5) becomes
(8)Di+1/2,j,kSi+1,j,kn+1−Si,j,kn+1Δx2+Di−1/2,j,kSi−1,j,kn+1−Si,j,kn+1Δx2+Di,j+1/2,kSi,j+1,kn+1−Si,j,kn+1Δy2+Di,j−1/2,kSi,j−1,kn+1−Si,j,kn+1Δy2+Di,j,k+1/2Si,j,k+1n+1−Si,j,kn+1Δz2+Di,j,k−1/2Si,j,k−1n+1−Si,j,kn+1Δz2=ϕSi,j,kn+1−Si,j,knΔt,where
(9)Di+1/2,j,k=Di,j,k+Di+1,j,k2,Di−1/2,j,k=Di,j,k+Di−1,j,k2,Di,j+1/2,k=Di,j,k+Di,j+1,k2,Di,j−1/2,k=Di,j,k+Di,j−1,k2,Di,j,k+1/2=Di,j,k+Di,j,k+12,Di,j,k−1/2=Di,j,k+Di,j,k−12.Figure 3
Five-point difference grid.Then, a one-dimensional mathematical model of countercurrent imbibition was established. As Figure4 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:
(10)∂2Sw∂x2=φD∂Sw∂t,0≤x≤L,Sw0,t=Swm,x=0,t>0,∂Sw∂x=0,x=L,t>0,Swx,0=Swi,0<x<L,t=0.Figure 4
1D model used for the numerical solution.For the 2-n-1 grid, the difference equation is(11)Di+1/2Si+1n+1−Sin+1Δx2+Di−1/2Si−1n+1−Sin+1Δx2=ϕSin+1−SinΔt,where
(12)Di+1/2=Di+Di+12,Di−1/2=Di+Di−12.For the first grid,(13)Sw=Swm,S=0.For thenth grid,
(14)2Di−1/2Si−1n+1−Sin+1Δx2=ϕSin+1−SinΔt.To ensure the convergence of the iteration, that is, the algorithm’s stability, according to the Von Neumann condition of the parabolic equation,(15)r=DϕΔtΔx2≤12.For a givenΔx, Δt needs to be small enough to satisfy Equation (15) in order to ensure the stability of the numerical calculation.Figure5 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.Figure 5
Relative permeability and capillary pressure curves.Table 1
Typical low permeability core properties and fluid parameters.
ParameterValueParameterValueLength of core (cm)4Permeability (10-3 μm2)10Surface tension (mN/m)47.4Porosity (f)0.2Irreducible water (f)0.3Water viscosity (mPa·s)1Residual oil (f)0.3Oil viscosity (mPa·s)2Kr endpoint of water krwe0.1Water Coreynw4Kr endpoint of oil kroe0.65Oil Coreyno2Figure 6
Water saturation distribution in respect to position and time from the results of the numerical method.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:(16)R=S¯w−SwiSwm−Swi=1/L∫0LSwx,tdx−SwiSwm−Swi,whereL 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.


### 2.3. 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 coefficientdwSw is replaced by a constant d, Equation (3) becomes
(17)∂2Sw∂x2=ϕd∂Sw∂t.


#### 2.3.1. One-Dimensional Model

A one-dimensional model was constructed (Figure7) 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 Swi. The oil-water contact surface is located at x=0 and x=L, where the water saturation is Swm. Thus, according to the hypothesis, the following mathematical model can be obtained:
(18)∂2Sw∂x2=ϕd∂Sw∂t,0<x<L,Sw0,t=Swm,x=0,t>0,SwL,t=Swm,x=L,t>0,Swx,0=Swi,0≤x≤L,t=0.Figure 7
One-dimensional water imbibition model (the blue arrows represent the water flow, and the red arrows represent the oil flow).By separating the variables, the distribution law of the water saturation through the core can be obtained as follows:(19)Swx,t=Swm−4πSwm−Swi∑n=0∞12n+1e−2n+12π2tDsin2n+1πxD,n=0,1,2,⋯,where tD=dt/ϕL2,xD=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:
(20)S¯wx,t=Swm−8π2Swm−Swi∑n=0∞12n+12e−2n+12π2tD.The oil recovery of water imbibition att is
(21)R=S¯w−Swi1−Swi=Swm−Swi1−Swi1−8π2∑n=0∞12n+12e−2n+12π2tD.The cumulative oil productionQt can be calculated using Equation (21), and the derivative with respect to time can be taken to obtain the oil production:
(22)qt=dQtdt=8AdLSwm−Swi∑n=0∞e−2n+12π2tD.


#### 2.3.2. Verification

The determination of the average capillary diffusivity coefficientd 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 Figure5 and Table 1 to calculate the oil recovery of this model, and the result is shown as the dotted line in Figure 8. Then, we determined the integral mean value of the capillary diffusion coefficient:
(23)DA=∫Swi1−SorDSwdSw1−Sor−Swi.Figure 8
Comparison of the oil recovery results obtained using the analytical and numerical methods.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 dwSw. 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.


#### 2.3.3. Three-Dimensional Model

As is shown in Figure9, 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. Lx, Ly, and Lz 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.Figure 9
Three-dimensional water imbibition model (the blue arrows represent the water flow, and the red arrows represent the oil flow).The water saturation distribution of the countercurrent imbibition in the three-dimensional model is(24)Swx,y,z,t=Swm−64π3Swm−Swi∑n=0∞∑m=0∞∑k=0∞12n+12m+12k+1×e−2n+1π/Lx2+2m+1π/Ly2+2k+1π/Lz2dt/ϕ×sin2n+1πxLxsin2m+1πyLysin2k+1πzLz,n,m,k=0,1,2,⋯.The oil recovery of the water imbibition att is
(25)R=Swm−Swi1−Swi1−512π6∑n=0∞∑m=0∞∑k=0∞12n+122m+122k+12×e−2n+1π/Lx2+2m+1π/Ly2+2k+1π/Lz2dt/ϕ.The oil production of the water imbibition is(26)qt=512LxLyLzdπ4Swm−Swi∑n=0∞∑m=0∞∑k=0∞12n+122m+122k+12×e−2n+1π/Lx2+2m+1π/Ly2+2k+1π/Lz2dt/ϕ×2n+1Lx2+2m+1Ly2+2k+1Lz2.


### 2.4. Imbibition Oil Productivity considering the Facture Network

To predict the imbibition capacity of the fracture network (Figure10) 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 10
Fracture network model of volumetric fracturing near the wellbore.Figure11 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 fx and fy, 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).Figure 11
Top view of the fracture network model of a volumetric fractured well.Table 2
Categories of matrix blocks.
Matrix block typeNumberTwo faces in contact with fractures (1 in Figure10)4Three faces in contact with fractures (2 in Figure10)2fx+fy−2Four faces in contact with fractures (3 in Figure10)fx−1fy−1Total matrix blockfx+1fy+1Based 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 att considering the volumetric fracturing is
(27)R=Swm−Swi1−Swi1−64π4∑n=0∞∑m=0∞4f12n+122m+12e−2n+1π/2a/fx+12+2n+1π/2b/fy+12dt/ϕ+2fx−1f12n+122m+12e−2n+1π/a/fx+12+2n+1π/2b/fy+12dt/ϕ+2fy−1f12n+122m+12e−2n+1π/2a/fx+12+2n+1π/b/fy+12dt/ϕ+fx−1fy−1f12n+122m+12e−2n+1π/a/fx+12+2n+1π/b/fy+12dt/ϕn,m=0,1,2,⋯.The oil production of the water imbibition considering the volumetric fracturing is(28)qt=64abhdπ2Swm−Swi∑m=0∞∑n=0∞4f12n+122m+12e−2n+1π/2a/fx+12+2m+1π/2b/fy+12dtϕ×2n+12a/fx+12+2m+12b/fy+12+2fx−1f12n+122m+12e−2n+1π/a/fx+12+2m+1π/2b/fy+12dt/ϕ×2n+1a/fx+12+2m+12b/fy+12+2fy−1f12n+122m+12e−2n+1π/2a/fx+12+2m+1π/b/fy+12dt/ϕ×2n+12a/fx+12+2m+1b/fy+12+fx−1fy−112n+122m+12e−2n+1π/a/fx+12+2m+1π/b/fy+12dt/ϕ×2n+1a/fx+12+2m+1b/fy+12.


## 2.1. Mathematical Model of Water Countercurrent Imbibition

McWhorter et al. [29] first proposed the one-dimensional incompressible flow equation for two-phase fluids considering capillary pressure:
(1)qtϕAfw′Sw∂Sw∂x−1ϕ∂∂xdwSw∂Sw∂x+∂Sw∂t=0,where Sw is the water saturation, fw′Sw is the derivative of the water cut, Ф is the porosity; qt is the sum of the water and oil flow rates (m3); A is the area of the contact surface of the matrix (m2); x is the distance in the x direction (m); and t is the imbibition time (s). dwSw is called the capillary diffusivity coefficient (m2/s), which is expressed as
(2)dwSw=−Kkrwkrokrwμo+kroμwdPcdSw,where K is the absolute permeability (10-3μm2); kro and krw are the relative permeabilities of oil and water, respectively; μw and μo are the viscosities of water and oil (mPa·s), respectively; and Pc 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 qt is equal to 0 in this case. The flow equation for countercurrent imbibition can be written as
(3)1ϕ∂∂xdwSw∂Sw∂x=∂Sw∂t.As can be seen from Equation (2), the capillary diffusivity coefficient is a function of the water saturation, and the relationship curve for dw and Sw is shown in Figure 2. dwSw is a bell-shaped function, and dwSwi=dw1−Sor=0, where Swi is the initial water saturation and Sor is the residual oil saturation of the core. Between these two endpoints, there is a certain saturation where dw has the highest value, dmax.Figure 2
Relationship curve of the capillary diffusivity coefficient and water saturation (the dotted line is the average capillary diffusivity coefficient’s value).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.


## 2.2. 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 usedD to represent dwSw, which is a variable, and Equation (3) is rewritten as
(4)∂∂xD∂S∂x=ϕ∂S∂t,where S=Swm−Sw/Swm−Swi, which is the normalized water saturation, and Swm=1−Sor is the highest water saturation of the core.By extending Equation (4) into three-dimensional space, we obtained
(5)∂∂xDx∂S∂x+∂∂yDy∂S∂y+∂∂zDz∂S∂z=ϕ∂S∂t,where Dx, Dy, and Dz 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 difference method, the first order central difference of the variables x, y, and z were determined:
(6)∂f∂x=fi+1−fi−12Δx,∂f∂y=fj+1−fj−12Δy,∂f∂z=fk+1−fk−12Δz.In addition, the first-order forward difference was determined for the variablet:
(7)∂f∂t=fn+1−fnΔt=δfΔt.The five-point difference grid is shown in Figure3, and Equation (5) becomes
(8)Di+1/2,j,kSi+1,j,kn+1−Si,j,kn+1Δx2+Di−1/2,j,kSi−1,j,kn+1−Si,j,kn+1Δx2+Di,j+1/2,kSi,j+1,kn+1−Si,j,kn+1Δy2+Di,j−1/2,kSi,j−1,kn+1−Si,j,kn+1Δy2+Di,j,k+1/2Si,j,k+1n+1−Si,j,kn+1Δz2+Di,j,k−1/2Si,j,k−1n+1−Si,j,kn+1Δz2=ϕSi,j,kn+1−Si,j,knΔt,where
(9)Di+1/2,j,k=Di,j,k+Di+1,j,k2,Di−1/2,j,k=Di,j,k+Di−1,j,k2,Di,j+1/2,k=Di,j,k+Di,j+1,k2,Di,j−1/2,k=Di,j,k+Di,j−1,k2,Di,j,k+1/2=Di,j,k+Di,j,k+12,Di,j,k−1/2=Di,j,k+Di,j,k−12.Figure 3
Five-point difference grid.Then, a one-dimensional mathematical model of countercurrent imbibition was established. As Figure4 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:
(10)∂2Sw∂x2=φD∂Sw∂t,0≤x≤L,Sw0,t=Swm,x=0,t>0,∂Sw∂x=0,x=L,t>0,Swx,0=Swi,0<x<L,t=0.Figure 4
1D model used for the numerical solution.For the 2-n-1 grid, the difference equation is(11)Di+1/2Si+1n+1−Sin+1Δx2+Di−1/2Si−1n+1−Sin+1Δx2=ϕSin+1−SinΔt,where
(12)Di+1/2=Di+Di+12,Di−1/2=Di+Di−12.For the first grid,(13)Sw=Swm,S=0.For thenth grid,
(14)2Di−1/2Si−1n+1−Sin+1Δx2=ϕSin+1−SinΔt.To ensure the convergence of the iteration, that is, the algorithm’s stability, according to the Von Neumann condition of the parabolic equation,(15)r=DϕΔtΔx2≤12.For a givenΔx, Δt needs to be small enough to satisfy Equation (15) in order to ensure the stability of the numerical calculation.Figure5 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.Figure 5
Relative permeability and capillary pressure curves.Table 1
Typical low permeability core properties and fluid parameters.
ParameterValueParameterValueLength of core (cm)4Permeability (10-3 μm2)10Surface tension (mN/m)47.4Porosity (f)0.2Irreducible water (f)0.3Water viscosity (mPa·s)1Residual oil (f)0.3Oil viscosity (mPa·s)2Kr endpoint of water krwe0.1Water Coreynw4Kr endpoint of oil kroe0.65Oil Coreyno2Figure 6
Water saturation distribution in respect to position and time from the results of the numerical method.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:(16)R=S¯w−SwiSwm−Swi=1/L∫0LSwx,tdx−SwiSwm−Swi,whereL 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.


## 2.3. 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 coefficientdwSw is replaced by a constant d, Equation (3) becomes
(17)∂2Sw∂x2=ϕd∂Sw∂t.


### 2.3.1. One-Dimensional Model

A one-dimensional model was constructed (Figure7) 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 Swi. The oil-water contact surface is located at x=0 and x=L, where the water saturation is Swm. Thus, according to the hypothesis, the following mathematical model can be obtained:
(18)∂2Sw∂x2=ϕd∂Sw∂t,0<x<L,Sw0,t=Swm,x=0,t>0,SwL,t=Swm,x=L,t>0,Swx,0=Swi,0≤x≤L,t=0.Figure 7
One-dimensional water imbibition model (the blue arrows represent the water flow, and the red arrows represent the oil flow).By separating the variables, the distribution law of the water saturation through the core can be obtained as follows:(19)Swx,t=Swm−4πSwm−Swi∑n=0∞12n+1e−2n+12π2tDsin2n+1πxD,n=0,1,2,⋯,where tD=dt/ϕL2,xD=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:
(20)S¯wx,t=Swm−8π2Swm−Swi∑n=0∞12n+12e−2n+12π2tD.The oil recovery of water imbibition att is
(21)R=S¯w−Swi1−Swi=Swm−Swi1−Swi1−8π2∑n=0∞12n+12e−2n+12π2tD.The cumulative oil productionQt can be calculated using Equation (21), and the derivative with respect to time can be taken to obtain the oil production:
(22)qt=dQtdt=8AdLSwm−Swi∑n=0∞e−2n+12π2tD.


### 2.3.2. Verification

The determination of the average capillary diffusivity coefficientd 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 Figure5 and Table 1 to calculate the oil recovery of this model, and the result is shown as the dotted line in Figure 8. Then, we determined the integral mean value of the capillary diffusion coefficient:
(23)DA=∫Swi1−SorDSwdSw1−Sor−Swi.Figure 8
Comparison of the oil recovery results obtained using the analytical and numerical methods.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 dwSw. 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.


### 2.3.3. Three-Dimensional Model

As is shown in Figure9, 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. Lx, Ly, and Lz 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.Figure 9
Three-dimensional water imbibition model (the blue arrows represent the water flow, and the red arrows represent the oil flow).The water saturation distribution of the countercurrent imbibition in the three-dimensional model is(24)Swx,y,z,t=Swm−64π3Swm−Swi∑n=0∞∑m=0∞∑k=0∞12n+12m+12k+1×e−2n+1π/Lx2+2m+1π/Ly2+2k+1π/Lz2dt/ϕ×sin2n+1πxLxsin2m+1πyLysin2k+1πzLz,n,m,k=0,1,2,⋯.The oil recovery of the water imbibition att is
(25)R=Swm−Swi1−Swi1−512π6∑n=0∞∑m=0∞∑k=0∞12n+122m+122k+12×e−2n+1π/Lx2+2m+1π/Ly2+2k+1π/Lz2dt/ϕ.The oil production of the water imbibition is(26)qt=512LxLyLzdπ4Swm−Swi∑n=0∞∑m=0∞∑k=0∞12n+122m+122k+12×e−2n+1π/Lx2+2m+1π/Ly2+2k+1π/Lz2dt/ϕ×2n+1Lx2+2m+1Ly2+2k+1Lz2.


## 2.3.1. One-Dimensional Model

A one-dimensional model was constructed (Figure7) 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 Swi. The oil-water contact surface is located at x=0 and x=L, where the water saturation is Swm. Thus, according to the hypothesis, the following mathematical model can be obtained:
(18)∂2Sw∂x2=ϕd∂Sw∂t,0<x<L,Sw0,t=Swm,x=0,t>0,SwL,t=Swm,x=L,t>0,Swx,0=Swi,0≤x≤L,t=0.Figure 7
One-dimensional water imbibition model (the blue arrows represent the water flow, and the red arrows represent the oil flow).By separating the variables, the distribution law of the water saturation through the core can be obtained as follows:(19)Swx,t=Swm−4πSwm−Swi∑n=0∞12n+1e−2n+12π2tDsin2n+1πxD,n=0,1,2,⋯,where tD=dt/ϕL2,xD=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:
(20)S¯wx,t=Swm−8π2Swm−Swi∑n=0∞12n+12e−2n+12π2tD.The oil recovery of water imbibition att is
(21)R=S¯w−Swi1−Swi=Swm−Swi1−Swi1−8π2∑n=0∞12n+12e−2n+12π2tD.The cumulative oil productionQt can be calculated using Equation (21), and the derivative with respect to time can be taken to obtain the oil production:
(22)qt=dQtdt=8AdLSwm−Swi∑n=0∞e−2n+12π2tD.


## 2.3.2. Verification

The determination of the average capillary diffusivity coefficientd 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 Figure5 and Table 1 to calculate the oil recovery of this model, and the result is shown as the dotted line in Figure 8. Then, we determined the integral mean value of the capillary diffusion coefficient:
(23)DA=∫Swi1−SorDSwdSw1−Sor−Swi.Figure 8
Comparison of the oil recovery results obtained using the analytical and numerical methods.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 dwSw. 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.


## 2.3.3. Three-Dimensional Model

As is shown in Figure9, 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. Lx, Ly, and Lz 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.Figure 9
Three-dimensional water imbibition model (the blue arrows represent the water flow, and the red arrows represent the oil flow).The water saturation distribution of the countercurrent imbibition in the three-dimensional model is(24)Swx,y,z,t=Swm−64π3Swm−Swi∑n=0∞∑m=0∞∑k=0∞12n+12m+12k+1×e−2n+1π/Lx2+2m+1π/Ly2+2k+1π/Lz2dt/ϕ×sin2n+1πxLxsin2m+1πyLysin2k+1πzLz,n,m,k=0,1,2,⋯.The oil recovery of the water imbibition att is
(25)R=Swm−Swi1−Swi1−512π6∑n=0∞∑m=0∞∑k=0∞12n+122m+122k+12×e−2n+1π/Lx2+2m+1π/Ly2+2k+1π/Lz2dt/ϕ.The oil production of the water imbibition is(26)qt=512LxLyLzdπ4Swm−Swi∑n=0∞∑m=0∞∑k=0∞12n+122m+122k+12×e−2n+1π/Lx2+2m+1π/Ly2+2k+1π/Lz2dt/ϕ×2n+1Lx2+2m+1Ly2+2k+1Lz2.


## 2.4. Imbibition Oil Productivity considering the Facture Network

To predict the imbibition capacity of the fracture network (Figure10) 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 10
Fracture network model of volumetric fracturing near the wellbore.Figure11 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 fx and fy, 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).Figure 11
Top view of the fracture network model of a volumetric fractured well.Table 2
Categories of matrix blocks.
Matrix block typeNumberTwo faces in contact with fractures (1 in Figure10)4Three faces in contact with fractures (2 in Figure10)2fx+fy−2Four faces in contact with fractures (3 in Figure10)fx−1fy−1Total matrix blockfx+1fy+1Based 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 att considering the volumetric fracturing is
(27)R=Swm−Swi1−Swi1−64π4∑n=0∞∑m=0∞4f12n+122m+12e−2n+1π/2a/fx+12+2n+1π/2b/fy+12dt/ϕ+2fx−1f12n+122m+12e−2n+1π/a/fx+12+2n+1π/2b/fy+12dt/ϕ+2fy−1f12n+122m+12e−2n+1π/2a/fx+12+2n+1π/b/fy+12dt/ϕ+fx−1fy−1f12n+122m+12e−2n+1π/a/fx+12+2n+1π/b/fy+12dt/ϕn,m=0,1,2,⋯.The oil production of the water imbibition considering the volumetric fracturing is(28)qt=64abhdπ2Swm−Swi∑m=0∞∑n=0∞4f12n+122m+12e−2n+1π/2a/fx+12+2m+1π/2b/fy+12dtϕ×2n+12a/fx+12+2m+12b/fy+12+2fx−1f12n+122m+12e−2n+1π/a/fx+12+2m+1π/2b/fy+12dt/ϕ×2n+1a/fx+12+2m+12b/fy+12+2fy−1f12n+122m+12e−2n+1π/2a/fx+12+2m+1π/b/fy+12dt/ϕ×2n+12a/fx+12+2m+1b/fy+12+fx−1fy−112n+122m+12e−2n+1π/a/fx+12+2m+1π/b/fy+12dt/ϕ×2n+1a/fx+12+2m+1b/fy+12.


## 3. 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 is3.9×10−3μm2. The oil viscosity under the reservoir conditions is 2.3 mPa·s and the water viscosity 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−8m2/s). The average initial 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.Figures12 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.Figure 12
Relationship between the fracture length and the oil recovery factor.Figure 13
Relationship between the fracture width and the oil recovery factor.Figures14 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 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 production and the better the productivity of the water imbibition.Figure 14
Relationship between the fracture length and oil production.Figure 15
Relationship between the fracture width and oil production.


## 4. 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 coefficientd 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.



---

*Source: 6693359-2021-02-28.xml*