Study of Single Phase Mass Transfer between Matrix and Fracture in Tight Oil Reservoirs

In tight fractured reservoirs, oil in matrices is mainly explored due to mass transfer mechanisms during the pressure depletion process. In the modeling of mass transfer in fractured reservoirs using the dual porosity concept, the shape factor is the most important parameter and should be described accurately. However, the current shape factors are not suited for tight oil reservoir simulation because the characteristics of tight oil reservoirs are not taken into account. In order to solve this problem, a new mass transfer function for tight fractured oil reservoirs is proposed by introducing a new time-related correction factor which could consider not only the existence of the boundary layer in nano-microscale throats in tight porous media but also the heterogeneous pressure distribution in matrix blocks. In addition, special contact relations between matrix and fracture are included. The correction factor presented in this study is verified using the experimental data and numerical simulation results. Data analysis results demonstrate that the lower and slower the pressure propagation velocity, the longer the duration time of unsteady flow compared to conventional reservoirs. Therefore, in the calculation of mass transfer flow in tight oil reservoirs, the unsteady flow between fracture and matrix cannot be ignored.


Introduction
At present, the oil production process of tight oil reservoirs is normally applied after the artificial fracturing procedure.During the pressure depletion process, oil in matrix blocks is produced by the mechanism of single phase mass transfer, normally.Similar to the conventional fractured oil reservoir, dual porosity is usually used to represent fractured reservoirs which was firstly proposed by Barenblatt et al. [1].Barenblatt et al. assumed that the flow from the matrix to the fracture is pseudo-flow; Warren and Root [2] proposed a pseudo-steady analytical radical solution of single phase mass transfer flow in fractured reservoirs and applied the analytical solution in well test analysis.The continuity equation describing the planar flow of a compressible fluid in a fracture is described as follows: On the assumption that the flow of the matrix follows Darcy's law, the rate expression of mass transfer is described by In Eq. ( 2), σ is the characteristic coefficient of the fractured matrix block and it is usually termed as the shape factor.Different scholars have great differences on the selection of shape factor values.In Warren and Root's model, the shape factor is defined as follows: In Eq. (3), n n = 1,2,3 represents the dimension of the matrix block.For the case of 1, 2, and 3 sets of fractures (1-D, 2-D, and 3-D situations), the values of σ are 12/L 2 , 32/L 2 , 60/ L 2 .Kazemi et al. [3] also have proposed the shape factor as follows in Eq. ( 4): In Eq. ( 4) L x , L y , and L z are the lengths of the matrix block in the x, y, and z directions, respectively.When L x = L y = L z = L, for the situation of 1, 2, and 3 sets of fractures, values of σ in different situations are 4/L 2 , 8/L 2 , and 12/L 2 .
Coats [4] has proposed different expressions of shape factor and compared with the fine grid numerical simulation.Simulation results reported by Bourbiaux shows that Coat's calculation results are more in coincidence with the reference solution, and Kazemi et al.'s results deviate from fine grid results.These mass transfer models mentioned above are proper for the conventional fractured oil reservoir in which the fluid flow in the matrix obeys Darcy's law and the mass transfer between the matrix and the fracture can be assumed as a steady or pseudo-steady flow.Hassanzadeh and Pooladi-Darvish [5] analyzed the effects of fracture boundary conditions on the matrixfracture transfer shape factor.Hassanzadeh and Pooladi [6] extend their previous analysis and use infinite-acting radial and linear dual-porosity models, where the boundary condition is chosen at the wellbore, as opposed to that at the matrix boundary.
Ranjbar and Hassanzadeh [7] used the matrix-fracture transfer shape factor for modeling the flow of a compressible fluid in dual-porosity media.Ranjbar et al. [8] had investigated the effect of the fracture pressure depletion regime on the shape factor for a single-phase flow of a compressible fluid.In the current study, a model for evaluation of the shape factor is derived using solutions of a nonlinear diffusivity equation subject to different pressure depletion regimes.A combination of the heat integral method, the method of moments, and Duhamel's theorem is used to solve this nonlinear equation.Ranjbar [9] analyzed one-dimensional matrix-fracture transfer in dual porosity systems with variable block size distribution.Ranjbar et al. [10] present a semianalytical solution for release of a single-phase liquid or gas from cylindrical and spherical matrix blocks with various block size distributions and different pressure depletion regimes in the fracture.However, at present, there is no mass transfer function for tight oil reservoirs.The main reason is that as for tight oil reservoirs, there are some particular characteristics compared to other reservoirs.These characteristics are due to the following: (1) a large amount of nanoscale or microscale throats present in tight porous media, with a boundary layer effect which leads to nonlinear flow behavior in matrix blocks [11,12], cannot be neglected; (2) the flow pattern has its own characteristics due to the contact relation between the matrix and the artificial fractures in tight oil reservoirs; (3) due to the low permeability of matrix blocks, the pressure sweep velocity is slow, the distribution of pressure in the matrix block is heterogeneous, and therefore, it is not appropriate to assume that the pressure distribution in the matrix is homogeneous and to use an average pressure at the center of the matrix to represent the pressure in the matrix.
In view of the abovementioned problems, a new mass transfer model is proposed.Firstly, in order to reflect the nonlinear flow behavior in matrix blocks, a model of tight formation permeability is used to replace the permeability in the traditional mass transfer model.Secondly, according to the contact relation between the matrix and the fracture, the new mass transfer model is divided into three categories which is the matrix-planar fracture model (1-D), the matrix-planar/naturally fracture model (2-D), and the matrixvolume fracture model (3-D).Then, based on the assumption of an unsteady flow between matrix blocks and fractures, a time-dependent correction factor is obtained to modify the traditional mass transfer flow model.At last, the proposed model is verified by comparing with experimental data or simulation results.

Mathematical Methods
2.1.Permeability Model of Tight Formation.The range of the tight oil formation throat radius distribution ranges from 20 nm to 1.2 μm, and the permeability is at the magnitude of 10 −1 mD [13].The existence of the boundary layer cannot be ignored in these nanoscale or microscale throats.The influence of the boundary layer on the throat is shown in Figure 1.The fluid flow in the throat of unconventional oil reservoirs reduces the effective flow due to the presence of the boundary layer.In some extreme cases, when the boundary layer thickness is equal to the original throat radius, the presence of the boundary layer can even cause all fluids in the throat to become immobile.According to Tian et al. [14], the boundary layer thickness can be quantitatively expressed by h is the thickness of the boundary layer, in μm, r is the radius of the throat in tight porous media, in μm, μ is viscosity, and ∇p is the pressure gradient, in MPa/m.

Geofluids
Due to the influence of the boundary layer, the radius of the original throat is larger than the effective flow radius.Therefore, the effective throat radius is used to characterize the flow in micro-or nanoscale throats.The effective throat radius is equal to the original throat radius minus the boundary layer thickness given by r eff = r i − h 6 The matrix core can be equivalent to a microcircular tube capillary bundle model.It is assumed that the model is composed of a circular tube with a continuous radius and that the flow of fluid in the matrix can be obtained according to the Poiseuille equation: For the pore throat distribution of a tight core, the Gauss probability distribution function (Eq.( 8)) can be used [15][16][17][18].Figure 2 shows the pore throat radius distribution and the cumulative pore throat distribution curve expressed by the Gauss function.Because of the different distribution of the pore throat in different cores, the pore throat distribution can be changed by changing the mean pore throat value and the standard deviation between pore throat and pore, thus representing the pore throat composition of different cores.For the distribution of pores and throats as normal distribution, the throat distribution can be described by the Gauss function, and the pore throat that cannot be fitted to the normal distribution can be represented by finding other functions.
Combining Eq. ( 7) and Eq. ( 8), the flow rate when the pores and throat distribution obey the Gauss function can be presented by Based on Darcy's law, the flow rate of porous media can be expressed as The matrix cross-sectional area can be expressed as The flow rate expression is obtained by substituting Eq. ( 11) into Eq.( 10): Eq. ( 12) and Eq. ( 9) are combined to obtain the permeability of the tight oil reservoir: 3 Geofluids Eq. ( 13) does not consider the effect of the boundary layer.Since the boundary layer effect cannot be neglected in the tight oil reservoir, the thickness of the boundary layer shall be subtracted from the original radius of the throats.Then, the permeability of the tight oil reservoir is obtained as 2.2.Mass Transfer Model.Because of the stress distribution, fracture network complexity, and the difference of fracturing technique, different forms of matrix fracture contact relations are formed in the reservoir after the fracturing procedure.Figure 3  The continuity equation of fracture can be expressed as The continuity equation of the matrix can be expressed as The mass transfer rate is a function of pressure difference, matrix permeability, and shape of matrix blocks, and it can be expressed as As shown in Figure 5, single phase mass transfer depends on the pressure difference between matrix blocks and fractures.The mass transfer rate can also be expressed by writing Darcy's law between the matrix and the fracture as given by q = Ak m μB The flow rate at the interface can be expressed as follows: Assuming a planar fracture, the area of interface is given by To calculate the mass transfer rate properly, the following assumptions are made: (1) The fracture pressure is chosen as the pressure at the interface p f due to the smaller size of the fracture relative to the matrix (2) The matrix pressure is selected as average pressure p m of the matrix block which is involved in the flow, and the position of average pressure p m is selected at the half width of the matrix that is involved in the flow (3) When 1-D single phase mass transfer occurred, the position change of the average pressure of the matrix is shown in Figure 5: Due to the extremely low permeability of the tight matrix block, mass transfer between matrix and fracture is slow.During every step of this process, only part of the matrix block is involved to contribute to the mass transfer (shown in Figure 5).As the mass transfer process proceeds, more and more matrix block regions are involved, and the average pressure point is selected at the center of the involved matrix block which contributes to the mass transfer, so the average pressure point (red point in Figure 5) continually changes its position as the mass transfer process proceeds.
4 Geofluids These assumptions are not consistent with the actual situation, but convenient for calculation.To alleviate this problem, a correction factor C f is introduced and then the mass transfer rate can be expressed as For a single phase slightly compressible fluid, the flow in the matrix can be written using subject to the initial condition and the boundary condition 5 Geofluids If M t is the cumulative mass transfer from matrix to fracture at time t, then M ∞ represents the total mass transfer when the time approaches to infinity; the ratio of M t and M ∞ can be given as (Crank, 1995) The ratio given by Eq. ( 26) can also be transformed into the ratio of increment of density: Fluid is assumed as slightly compressible, thus Combining Eq. ( 17)-Eq.( 28), an analytical solution is obtained as follows: Because the mass transfer rate is equal to the accumulate flow per unit volume of matrix, The partial derivative of pressure P m to time t is as follows: Substituting Eq. (31) into Eq.( 22), In Eq. (33), A 1 = 2L 2 , V is the volume of the matrix involved in the mass transfer process and the correction coefficient can be expressed as follows: Substituting Eq. ( 32) and the pressure difference of the average pressure of the matrix and the pressure of the fracture (P m − P f ) into Eq.( 34), the final expression of correction coefficient C f is obtained as where t D = k m /ϕ m μc t L 2 t When the value of dimensionless time t D is large enough, the correction factor converged into the shape factor of the pseudo-steady mass transfer flow: For the correction coefficient C f of pseudo-steady mass transfer flow, Lim and Aziz use the constant shape factor σ to replace C f ; the mass transfer rate can be expressed as follows: where V = LA 0 is the volume of the matrix and A = 2A 0 , The correction coefficient of the pseudo-steady mass transfer coefficient is then  6 Geofluids pressure in the matrix during the process of mass transfer is also shown in Figure 7, where after the simplification the contact relation becomes a circular contact between matrix and fracture.Under the assumption of 2-D circular contact relation, the pressure diffusion equation in the matrix block can be written in the following form:

∂P ∂r 40
The boundary and initial conditions are Then, for the analytical solution of two-dimensional radical flows, after handing, the pressure difference relation can be written as In Eq. ( 43), J 0 is the Bessel function of the first kind of order zero, and α n is the root of the Bessel function.
After substituting Eq. (44) into Eq.( 45), the mass transfer rate of Eq. ( 46) can be obtained as given by Figure 7 shows that the position of the average pressure in 2-D situation changes as the mass transfer process keeps proceeding.The theory is the same as in Figure 5.
In Eq. ( 47),   Under the assumption of 3-D contact, the pressure diffusion equation in the matrix block can be written in the following form:

∂P ∂r 49
The boundary condition and initial condition are The analytical solution of the 3-D pressure diffusion equation can be written as follows [19]: The correction factor of the mass transfer equation in 3-D is In Eq. (52),

Geofluids
Eq. ( 53) shows the correction factor in the circumstance of 3-D and in Eq. ( 53)

Model Validation
3.1.Validation of Tight Oil Permeability Model.For the validation of the tight oil permeability model, the accuracy of the model is verified by comparing with the experimental results obtained by Wang et al. [13,20].In his experiment, throat size distribution and core permeability are obtained, respectively, from mercury injection and the core displacement experiment.
Main parameters and throat size distribution used in the core experiment are listed in Tables 1 and 2. Based on the distribution of throat size and some other related reservoir parameters, the effective permeability of the tight oil reservoir can be calculated from the proposed permeability model.A comparison of the calculated and experimental results is shown in Figure 9.The error line gives the difference between the experimental and measured values.The results shown in Figure 9 demonstrate an acceptable agreement between experimental and calculated permeabilities which verifies the tight oil reservoir permeability model.

Validation of Tight Oil Mass Transfer Model.
For the matrix-planar fracture and matrix-planar/naturally fracture model, their corresponding correction factors are compared with the constant shape factor which is previously proposed by Lim and Aziz [21].The main parameters used in calculation are listed in Table 3, and the change of the mass transfer flow correction factor with time is shown in Figure 10.The results are shown in Figure 11 and demonstrate that the tight oil mass transfer correction factor decreases as time increases and finally tends to a stable value similar to the traditional shape factor [21].
For formation with different permeabilities, the correction factors of mass transfer flow are different (Figure 10).For the tight oil reservoir, since the permeability of the reservoir is low, the pressure propagation is low and the time at which the correction factor reaches steadiness is obviously longer than that of high permeability reservoirs.Therefore, when calculating the mass transfer flow in tight oil reservoirs, the time-related correction factor should be considered.For high-permeability reservoirs, due to highpressure diffusivity, the time of the correction coefficient reaches steadiness much faster than that of tight oil reservoirs.Therefore, the time-related correction coefficient has great influence on the mass transfer flow in highpermeability reservoirs.
For the 3-D model, the mass transfer rate and cumulative flow are compared with numerical simulation results; the traditional mass transfer model with a constant shape factor is previously proposed by Lim and Aziz [21].The main parameters used in calculation are also listed in Table 3, except the permeability of the matrix which is 0.1 mD.
A comparison of the mass transfer rate is shown in Figure 11(a), and the results show the mass transfer rate of the tight oil reservoir which agrees well with numerical simulation results except the very early stage of the mass transfer process.This could be caused by the assumption of radical flow which may not coincide with the real situation at the early stage.Besides, in Figure 11(a) the mass transfer rate with a constant shape factor is obviously lower than the one with the time-related correction factor at the early and middle stages.However, during the late stage, as the time-related correction factor changes, its value gradually approaches to the value of the constant shape factor, then the mass transfer rate with the time-related correction factor turns the same with the mass transfer rate with the constant shape factor.Due to the better fitting of numerical simulation results, the proposed correction factor C f is proven to be more suitable for the tight oil matrix block.On the other hand, the corresponding cumulative mass transfer is shown in Figure 11(b) and compared with the constant shape factor.The cumulative mass transfer with the time-dependent correction factor is much close to the numerical simulation results due to the higher mass transfer rate at the early and middle stages of the mass transfer process.Overall, even though there is little discrepancy, the result of the new model    Mean pore radius, μm µ:

Conclusion
Oil viscosity, mPa•s ▽p: Pressure drop, MPa/m r eff : Effective throat radius, μm r i : Original throat radius, μm f(r): The micro throat distribution function N: The total number of microtubes A: The cross-sectional area of the core, μm 2 φ m : Porosity of matrix, % φ f : Porosity of fracture, % Q: Total flow rate of capillary bundle model, cm 3 /day L: Length of capillary bundle model, cm k m : Matrix permeability of tight oil reservoir, mD k f : Fracture permeability, mD p f : Average pressure in fracture, mD p m : Average pressure in matrix block, mD q m→f : Mass transfer flow rate, cm 3 /day C t : Total compressibility Q: Mass transfer rate, cm 3 /day V f : Volume of fracture, cm 3 M ∞ : Total mass transfer quantity when the time approaches to the ultimate time M t : Accumulate mass transfer quantity t at time t ρ i : Oil density at initial time, g/cm 3 ρ f : Oil density in fracture at time t, g/cm 3 ρ m : Average oil density in matrix, g/cm 3 C f : Correction factor, dimensionless.

Figure 1 :Figure 2 :
Figure 1: The influence of the boundary layer on the flow of tight throat.
(a) shows a matrix-planar fracture contact model in 1-D, Figure 3(b) shows the matrix-planar/naturally fractured model in 2-D, and Figure 3(c) shows the matrix-complex fracture in 3-D.For a different contact model of matrix and fractures, mass transfer will happen when the original balanced pressure distribution changes.Analytical solutions of the mass transfer rate of different contact relation models are characterized in the following section.

2. 2 . 1 .
Matrix-Planar Fracture Mass Transfer Model.The contact relation of the matrix and fracture in 1-D is shown in Figure4.Under the pressure difference, fluid in the matrix flows linearly into the fracture.During the process of linear flow, the pressure in the matrix changes continually.Thus, in the calculation process, pressure change in the matrix should be considered to avoid introduction of a significant error.

Figure 6 FractureFigure 5 :
Figure 5: Schematic for the change of matrix average pressure.

Figure 7 :Figure 6 :
Figure 7: Schematic for the change of average pressure.

Figure 8 :
Figure 8: Schematic for contact between matrix and volume fractures.
Figure 8   shows the schematic of contact between matrix blocks and volume fracture in 3-D.Under this circumstance, matrix blocks are surrounded by fractures, and fluid flow from matrix blocks to fracture can be simplified as spherical flow.

휇m 2 )Figure 9 :
Figure 9: Comparison between experimental data and calculated data from the model.

Figure 10 :
Figure 10: The change of the mass transfer flow correction factor with time.

9
Geofluidsis much closer to the real situation of the mass transfer in the tight fractured oil reservoir.

( 1 ) 2 )( 4 )
A tight matrix permeability model is established by considering both the boundary layer and the throat distribution of tight oil reservoirs.A nonlinear flow in the matrix can be described through the new matrix permeability model (Three kinds of models have been proposed to reflect the contact relation of matrix and fracture based on the actual complex fracture network distribution during the fracturing process: matrix-planar fracture (1-D) model, matrix-planar/naturally fracture (2-D) model, and matrix-volume fracture (3-D) model (3) The correction factor has considered two main characteristics in tight fractured oil reservoirs: the first is the change of pressure distribution in the matrix during the mass transfer process and the second is the existence of the boundary layer in tight porous media.Data analysis results show that despite little discrepancy at the early stage, the newly proposed correction factor made a much more accurate calculation result of the mass transfer rate compared to the traditional mass transfer function with a constant shape factor For fractured tight oil reservoirs, the pressure propagation velocity in the matrix is low, and the duration time of the unsteady mass transfer flow between matrix fractures is very long, so the unsteady mass transfer flow plays an important role in the exploration of tight oil reservoirs Nomenclature H: Thickness of boundary layer, μm R: Radius of throat, μm υ:

Table 2 :
Data of throat distributions of Yanchang Formation.

Table 3 :
The main parameters of the model.