A Semianalytical Two-Phase Imbibition Model in Dual-Porosity and Dual-Permeability Reservoirs

The pressure transient behavior of water injection well has been extensively investigated under single-phase flow conditions. However, when water is injected into formation, there are saturation gradients within the water flooded area. Additionally, water imbibition is essentially important for oil displacement in dual-porosity and dual-permeability (DPDP) reservoirs. In this work, a novel semianalytical two-phase flow DPDP well test model considering both saturation gradient and water imbibition has been developed. The model was solved by the Laplace transform finite difference method. Type curves were generated, and flow regimes were identified by the model. The model features and effect of parameters were analyzed. Results show that water imbibition reduces the advancing speed of water drive front in the fracture system and slows down the water cut raising rate and the expansion speed of the two-phase zone in the fracture system. Therefore, the fluid exchange between the fracture and matrix systems becomes more sufficient and more oil will be recovered from the DPDP reservoir. The shape of pressure curves is similar for the single-phase and two-phase flow DPDP model, but the position of the proposed model is above the curves of the singlephase model. Shape factor mainly influences the interporosity period of the pressure derivatives. Water imbibition has a major effect on the whole system radial flow period of the curves. The findings of this study can help for better understanding of the oil/water two-phase flow pressure transient behavior in DPDP reservoirs considering saturation gradients and water imbibition.


Introduction
Fluid flows in natural fractured reservoirs are usually described by a dual-porosity media. The dual-porosity systems are two overlapping continua with uniform and orthogonal fractures and matrix. The fractures provide the conductive pathway for fluid flow, whereas the matrix is treated as a source term to fractures [1,2]. If both the fracture and matrix have continuity and the matrix conductivity is significant, the model is considered as the dual-permeability model [3]. Liu et al. first proposed the exact solution of a single-phase DPDP model with consideration of wellbore storage and skin effect [4,5]. After that, several analytical solutions [6][7][8][9][10][11] for the DPDP single-phase well test model were derived with different well and formation conditions. Guo et al. [12] proposed a triple-porosity and dualpermeability pressure transient model for a horizontal well in fractured-vuggy carbonate reservoirs. The vugs are assumed to be dispersed throughout the reservoir, and the media connected to the horizontal wellbore are treated as a DPDP system. Liu et al. [13] presented an analytical trilinear flow model for a multiple-fractured horizontal shale gas well. The model divided the reservoir into three regions: hydraulic-fracture region, microfracture network, and pure-matrix region. The microfracture is treated as a dual-porosity media. However, the permeability in the matrix is neglected in the dual-porosity system. Nie et al. [14] developed a semianalytical DPDP model for a horizontal well and obtained the type curves by numerical simulation method. It indicates that the single-permeability modeling ignores the direct fluid supply from the matrix to the wellbore, resulting in an obvious difference between the dual-permeability and single-permeability models. And there are also some numerical DPDP models used for the simulation of hydraulic fractures and natural fractured reservoirs [15][16][17][18][19][20].
The methods for well test analysis of water injection well mostly still use the single-phase well test method. However, when water is injected into the formation, the flow becomes oil/water two-phase flow during water injection. Therefore, it is necessary to investigate the two-phase flow pressure behavior in the DPDP model for water injection well [21]. The reservoir is divided into two or three fluid banks (the flooded region, the displaced fluid bank, and the unflooded region) during water flooding [22][23][24][25][26]. A saturation gradient is formed after water flooding due to the differences in oil and water properties. The location of the water flood front is determined based on the Buckley-Leverett (1942) frontal advance equations [27,28]. The water saturation distribution can be estimated by discretizing the fluid bank into a series of banks [29][30][31].
Based on Aronofsky's classical imbibition exponential equation [32], De Swaan proposed the equations for water imbibition rate in a fracture surrounded by a matrix in a convolution form [33]. By combining the Buckley-Leverett theory with the imbibition model, the water saturation distribution considering imbibition can be estimated for water injection well in fractured reservoirs [34,35].
According to the quasistationary method [36], water saturation can be regarded as constants when solving the diffusivity equation for pressure. Thus, the water saturation is decoupled from the pressure, and the pressure can be solved by the Laplace transform finite difference (LTFD) method in the Laplace domain. This method is semianalytical in time. Hence, it eliminates the stability and the convergence issues resulting from time discretization [37,38].
The models for pressure transient analysis in the DPDP reservoir mostly still use a single-phase method. The twophase flow model considering the water saturation change in the literatures usually divided the reservoir into a series of banks, and the model is solved by the multicomposite method [39], see Figure 1. The quasistationary method and LTFD method are rarely applied in the solution of a twophase flow model. In this work, we proposed a two-phase flow model for a water injection well in the DPDP reservoir. First, we developed the water saturation model considering water imbibition based on Kazemi's linear imbibition model [34]; then, the water saturation model is coupled with the press transient model according to the quasistationary method. The semianalytical model is solved by the LTFD method. Model features and effects of the parameters on pressure behavior of the DPDP model were also discussed in detail. According to this study, it is necessary to consider the two-phase flow pressure behavior for water injection in natural fractured reservoirs, which improves the accuracy of well test interpretation results.

Model Description
Based on the quasistationary method, water saturation is decoupled from the diffusion equation. Therefore, we first derived the water saturation model for water injection well in the DPDP reservoir. Besides, water imbibition is also considered in the saturation model. Then, the pressure transient DPDP model for water injection well considering saturation gradients was developed and solved by the LTFD method. For the two-phase proposed model, the parameters of total compressibility, interporosity flow coefficient, and total mobility change with water saturation, which is different from the conventional two-phase well test method based on the P-M theory [40,41].

Water Saturation Model
2.1.1. Physical Model. The natural fractured reservoir is assumed to be a horizontal radial DPDP reservoir with an upper and lower sealed boundary. The reservoir has a uniform initial pressure and constant temperature. The matrix system has storage capacity and conductivities, and fluid flow exists in both the matrix and fracture systems. The fluids flow through the fracture and matrix are slightly compressible and obey Darcy's law. Water imbibition and the mass exchange between matrix and fracture are considered. When water is injected into the DPDP formation, it is divided radially into two regions by a moving interface. As seen in Figure 2, r w , r f , and r e indicate the radius of the wellbore, the position of water drive front, and the radius of the reservoir outer boundary, respectively. Region 1 is the water flooded region, and region 2 is the unflooded region.  2 Geofluids cumulative oil recovery by imbibition from a piece of rock surrounded by water can be written in an exponential form [32]: where Imbibition in a fracture surrounded by matrix was also investigated by De Swaan, and he calculated the rate of imbibition in a convolution form based on a water volume balance in a fracture [33]. Kazemi et al. deduced a liner flow imbibition model that accounts for saturation changes in fracture and matrix [34]. In this section, we extend the formula to a radial flow system for the DPDP model. The water saturation change in the fracture can be expressed as And the water saturation in the matrix can be written as Initial conditions Boundary conditions Introducing Laplace transform, the water saturation in the Laplace domain can be written as where z 1 is the Laplace variable, The water saturation in fracture and matrix in the real domain can be obtained by inversion of Laplace transform: The water saturation distribution in the matrix and fracture systems can be calculated according to the analytical solution of the water saturation model. The basic parameters for the saturation model calculation are listed in Table 1. Figures 3 and 4 are the water saturation curves in the fracture and matrix systems changing with distance, when the imbibition rate coefficients are 0.01 and 1, respectively.  Figure 3 illustrates the water saturation distribution in the fracture and matrix systems varying with distance when the imbibition rate coefficient is 0.01 per day. The solid lines marked with the "×" symbol are water saturation curves in the fracture system, and the unmarked solid lines are water saturation curves in the matrix system. The red, green, light blue, pink, and blue colors indicate the water injection time is 5, 50, 100, 200, and 300 days, respectively. From the figure, when the imbibition rate coefficient is small, the oil/water two-phase zone becomes larger with the increase of injection time. The advancing speed of water drive front in the fracture system is faster than that in the matrix system. The reason is

Geofluids
that only a small amount of water flows into the matrix by water imbibition due to the small value of the imbibition rate coefficient. The fluid exchange is not enough between the fracture and matrix systems. Thus, only part of the crude oil in the matrix system is displaced into the fracture system and the two-phase zone will be enlarged. Figure 4 illustrates the water saturation distribution in the fracture and matrix systems varying with distance when the imbibition rate coefficient is 1 per day. The solid lines marked with the "×" symbol are water saturation curves in the fracture system, and the unmarked solid lines are water saturation curves in the matrix system. The red, green, light blue, pink, and blue colors indicate the water injection time is 5, 50, 100, 200, and 300 days, respectively. By comparison of the curves, it is found that when the imbibition rate coefficient becomes larger, more water is exchanged into the fracture system and the advancing speed of water drive front in the fracture system is basically the same with that in the matrix system. Water imbibition prohibits water coning and slows down the water cut raising rate in the fracture system. Therefore, water imbibition reduces the expansion speed of the two-phase zone and improves the oil recovery in the natural fractured reservoir.

Pressure Transient Model.
The pressure transient DPDP model takes into account the fluid flow in the matrix. The interporosity flow between the matrix and fracture systems is assumed to be pseudosteady. According to the quasistationary method, the saturation change in time can be considered invariant when the diffusivity equation is solved for pressure. Thus, the radial diffusion equations of the DPDP model are transformed into the Laplace domain and solved by the finite difference method. The diffusion equations for the two-phase flow DPDP model are presented in Appendix A.
A logarithmic transform, z = ln ðr D Þ, is used to change the unequal radial grid into equal distant radial grid. Taking the Laplace transform with the quasistationary assumption, the diffusion equations for the two-phase flow DPDP model can be written as The initial conditions The inner boundary conditions considering the skin effect The inner boundary conditions considering the wellborestorage effect The constant pressure outer boundary conditions The no-flow outer boundary conditions By using the Laplace transform finite difference approximation, the formation is divided into n grid blocks in radial direction according to the point-centered logarithmic gridding method. The diffusion equations and boundary conditions can be approximated by the LTFD method, and the detailed derivations can be seen in Appendix B. The pressure difference equations can be written in a matrix form as where M is a large sparse matrix, x * is the unknown pressure By computing the coefficient relation matrix, the bottomhole pressure of the water injection well in the Laplace domain is obtained. Based on the Stehfest algorithm [42], 6 Geofluids the dimensionless wellbore pressure and pressure derivative of the water injection well for the proposed model can be solved. The basic parameters and the relative permeability data for the calculation of pressure curves are listed in Table 2 and Figure 5.

Results and Discussion
3.1. Model Validation. The model solutions are validated by a field well test in a DPDP reservoir. The tested well is a water injection well with a constant water injection rate. We obtain the desired interpretation results of the well based on the proposed model by adjusting the parameters repeatedly and doing fitting experiments continuously. The basic parameters for the water injection well are listed in Table 3. Figure 6 shows the pressure behavior of the injection period for water injection well in the DPDP reservoir. As shown, there is a good agreement between the proposed model and the field test. It indicates that the model in this work is capable of providing a reliable basis for field dynamic analysis.

Model Feature Analysis.
The pressure and the pressure derivative curves of the two-phase flow DPDP model are plotted for model feature analysis. The flow regimes of type curves are identified. The effect of the parameters on the type curves is also analyzed in this section. Figure 7 illustrates the pressure transient behavior of the proposed model without considering water imbibition. "P wD " and "dP wD " in the figure represent the dimensionless pressure and pressure derivative of the bottom-hole, respectively. "t D " is the natural logarithm  of dimensionless time. According to the feature of the curves, the formation flow can be divided into four stages. The first stage is the early flow period due to the wellbore storage effect and skin effect. The second stage is the interporosity flow period caused by the pseudosteady flow between the fracture and the matrix. The third stage is the transition flow period between the second and fourth stages. The fourth stage is the whole system, including the fracture and matrix systems, flow period. The horizontal line is the characteristics for the whole system radial flow period.

3.2.2.
Comparison between the Single-Phase and the Proposed Two-Phase Models. Figure 8 shows the pressure curves of the single-phase DPDP model and the proposed model without considering water imbibition. The basic parameters for the single-phase model are the same as the proposed model. The single-phase model is also solved by the LTFD method. From the figure, the shape of the curves for the proposed model and single-phase model is similar, but the position of the proposed model moves upward from the position of the single-phase model. The reason is that the two-phase fluid properties, such as fluid viscosity, density, relative permeability, and total mobility, are lower than that of the single phase. Therefore, the flow capacity decreases and the pressure curves go upward when it becomes a two-phase flow.    Figure 9: Influence of different shape factors on pressure transient behavior. 8 Geofluids Figure 9 describes the effect of shape factor on pressure curves of the proposed model. The red, green, pink, and blue colors indicate the shape factor values of 10 -6 , 10 -5 , 10 -4-, and10 -3 , respectively. From the figure, it is concluded that the shape factor has major influence on the interporosity period of the pressure derivatives. With the increase of the shape factor, the duration of the interporosity period becomes smaller. When the value of the shape factor gets greater, the interporosity flow between the fracture and the matrix becomes easier. Thus, the interporosity period becomes longer as the shape factor decreases.

Effect of Imbibition
Rate Coefficient. Figure 10 shows the pressure transient behavior of the DPDP model when the imbibition rate coefficient values are 0, 0.001, 0.01, and 0.1, respectively. Water imbibition is characterized by the imbibition rate coefficient, which is a constant equal to the reciprocal of the time necessary to produce 63 percent of the oil recoverable from the rock [32]. It represents    9 Geofluids the capacity for the water to displace oil from the matrix with imbibition. The capacity becomes greater when the imbibition rate coefficient gets larger. From the figure, water imbibition mainly influences the whole system radial flow period of the pressure curves. The pressure derivative curve is a straight line at the whole system radial flow period when the imbibition rate coefficient is zero. As the imbibition rate coefficient increases, the pressure derivative curves rise up to a greater degree at the end of the radial flow period.

Effect of Boundary
Conditions. The influences of different boundary conditions on the pressure curves are illustrated in Figure 11. The red and blue curves represent the constant pressure and no-flow outer boundary conditions, respectively. By comparison of the curves, it is found that the pressure and pressure derivative curves rise up to a straight line due to the no-flow outer boundary. And the slope of the straight line is 1. The pressure derivative curve drops downward under constant pressure outer boundary.

Summary and Conclusions
In this work, a semianalytical two-phase flow model was proposed to analyze the flow behavior of the water injection well in DPDP reservoirs. The saturation gradients and water imbibition were included in the model, which was different from the conventional two-phase well test model. Based on the investigation, the following conclusions are drawn: (1) The possible flow regimes for the two-phase DPDP model include (1) early wellbore storage and skin effect period, (2) interporosity flow period, (3) transitional flow period, and (4) the whole system radial flow period (2) Water imbibition effectively controls the water breakthrough in the fracture system. With the increase of the imbibition rate coefficient, more water is imbibed from the fracture into the matrix. The oil in the matrix system is recovered mainly by water imbibition, while water flooding is the main displacement force in the fracture system (3) The shape of the two-phase pressure curves is similar to the single-phase curves, while the position of the two-phase curves is above the single-phase pressure curves. The pressure and pressure derivative curves rise up to a straight line with a slope of 1 due to noflow outer boundary and drop downward under constant pressure outer boundary (4) The shape factor mainly has an effect on the duration of the interporosity period. With the decrease of shape factor, the dip becomes flatter and the interporosity period gets longer. Water imbibition has major influence on the whole system radial flow period. With the increase of the imbibition rate coefficient, the pressure derivative curves rise up to a greater degree at the end of the radial flow period