Production Dynamic Characteristic of Fractured Wells in Multilayer Reservoirs Considering the Effect of Non-Uniform Flux Distribution

Hydraulic fracturing stimulation, which improves matrix permeability and reduces production costs, has been extensively used in the exploitation of multilayer reservoirs. However, little research on the production dynamic characteristics of vertically fractured wells in stratified reservoirs has been done in the literature. The influence of flux variation along the fracture on the pressure transient behavior has been ignored in these previous works. Therefore, this paper introduces a novel semi-analytical model for fractured wells in multilayer reservoirs, in which the finite difference method is used to characterize fluid flow in the fracture and the Green’s function method is used to characterize fluid flow in the matrix. With the aid of the model, the production dynamic characteristics of fractured wells in multilayer reservoirs can be readily investigated. In addition, based on the assumption of nonuniform flux distribution along the fracture, we successfully recognize four flow regimes occurring in the pressure drop and pressure derivative curves. Following that, the influences of several parameters on the pressure dynamics and layered flux contribution are studied. The calculation results indicate that a larger storability ratio, as well as a larger permeability ratio, can increase the values of the pressure drop and the pressure derivative; the greater the fracture height, the greater the fluid flow into each layer of the fracture. During the production of this model, increasing the fracture conductivity can reduce the pressure drop and pressure derivative, which means lower flow resistance in the fracture.With the growing dependence on fossil energy, deep and ultra-deep areas have steadily evolved into the next main potentials of resource exploration and development. Recently, China has consistently discovered a huge number of deep-layer reservoirs, such as the Puguang, Tahe, Shunbei, and Anyue oilfields, showing great resource potentialities and considerable economic benefits [1]. Considering the influence of the complex sedimentary environment, most deep reservoirs are composed of several layers with different stratigraphic characteristics. Commingling production is commonly adopted to increase producing profit for stratified reservoirs in oil and gas fields [2]. Given this, extensive literature was related to the pressure dynamics of a vertical well in stratified reservoirs [3-6]. Rahman and Mattar [7] derived a new analytical solution for the commingled-layered reservoir with unequal initial pressures in the Laplace domain. Onwunyili and Onyekonwu [8] developed a coupled model that can more accurately simulate the commingled production behavior of multilayer reservoirs. Shi et al. [9] investigated the impact of the vertical inhomogeneous closed boundary radii on pressure transient behaviors of the multilayered commingled reservoir. These previous researches give us a basic understanding of the production dynamic characteristics for the un-fractured stratified reservoirs.Nevertheless, in consideration of the large burial depth easily resulting in low permeability reservoirs, hydraulic fracturing stimulation has always been applied to improve the reservoir conductivity in the oilfield. However, the research on hydraulically fractured wells in multilayer reservoirs is short of detailed analyses. Several scholars have proposed analytical solutions related to fractured wells in multi-layer reservoirs [10, 11]. Bennett et al. [12] first presented analytical approximations and numerical solutions of vertically fractured wells without formation crossflow and correlated the multilayered solutions with the single-layer solutions. Osman [13] proposed a mathematical method to investigate the pressure responses of fracturing wells in stratified reservoirs and obtained theoretical curves of the model. Wang et al. [14] developed an analytical model of the vertically fractured well with infinite conductivity in boxed reservoirs and carried out the pressure transient analysis. Li et al. [15] presented a model to study the hydraulically fractured well with varying fracture conductivity and successfully investigated the influence of the fracture extension, fracture conductivity, storability factor, and transmissibility factor on pressure. However, the production dynamic characteristics of fractured wells in multilayer reservoirs is a quite complicated problem that needs to be considered in many aspects. In these analytical solutions, they oversimplify the case of flux distribution along the fracture by assuming the fluid flows uniformly into the fracture. For better matching the actual production field, it is necessary to investigate the production dynamic characteristics of a fractured well in layered reservoirs that consider the nonuniform flux distribution along the fracture.At present, this research proposes a semi-analytical method for studying the production dynamic characteristics of fractured wells in stratified reservoirs. With the aid of this model, finite fracture conductivity and flux variations along the fracture are considered in this paper. Subsequently, the accuracy of the method is authenticated by comparing the calculation results of commercial software (Eclipse). Ultimately, the influences of crucial parameters on the flux contribution of layers and pressure dynamics are discussed, including the wellbore storage factor, storability factor, permeability ratio, fracture height, fracture conductivity, and permeability anisotropy.To better analyze the production dynamic characteristics of hydraulically fractured wells in multilayer reservoirs, the physical model and mathematical model are introduced in this section. As shown in Figure 1, there is a vertical well located in the middle of the fracture in a multilayer reservoir. The vertical fracture cuts through both layers of the stratified reservoir. A single-phase fluid enters the wellbore only through the vertical fracture at a constant flux rate. Other assumptions are given below to construct the proposed semi-analytical model:The permeability, porosity, and production contribution for every layer are different;The matrix has homogeneous and isotropic characteristics;The upper and lower boundaries of the reservoir are impermeable;Each layer is separated by an interlayer, so the impact of crossflow should be ignored;The nonuniform flux distribution along the x-axis and z-axis is considered;The fluid properties are assumed to be slight compressibility and constant viscosity;The gravitational effect and capillary force are disregarded;The effect of the wellbore storage coefficient is considered.The permeability, porosity, and production contribution for every layer are different;The matrix has homogeneous and isotropic characteristics;The upper and lower boundaries of the reservoir are impermeable;Each layer is separated by an interlayer, so the impact of crossflow should be ignored;The nonuniform flux distribution along the x-axis and z-axis is considered;The fluid properties are assumed to be slight compressibility and constant viscosity;The gravitational effect and capillary force are disregarded;The effect of the wellbore storage coefficient is considered.Considering the fracture width is set to be much smaller compared to its length and height, The flow behavior occurring in the fracture can be simplified as two-dimensional flow. As shown in Figure 2, we discretize the fracture into Nf (Nf = Ni × Nj) panels along the x-direction and the z-direction. Note that the dimension of every panel is uniform, namely Δx × Δz. Furthermore, the flow behavior around a single fracture element is visualized in Figure 2. For the fracture system, the mathematic model that is not connected with the wellbore should be stated as follows [16]:where pf is the pressure of the fracture-panel, μ is the oil viscosity, qm-f is the flow from the matrix into the fracture-panel, β is a unit-conversion factor with a value of 0.0853, Δx and Δz are sizes of the fracture-panel w is the fracture-panel width, kf is fracture permeability, ϕf is fracture porosity, ctf is fracture total compressibility and t is time. With the aid of the dimensionless variables defined in Table 1, the dimensionless format of equation (1) can be obtained. Subsequently, applying the finite difference method to equation (1). For the nth time step, this equation can be rewritten as:where i and j is the position number of the fracture element. Here, the flow equation of the single fracture-panel (i, j) has been given. Applying equation (2) to all fracture elements, we can obtain Nf flow equations in the fracture system (see Supplementary Material 1).According to the assumption that the reservoir is separated by the interlayer, we know that the pressure and flux of each layer are discontinuous. Therefore, their pressure responses can be calculated respectively. Here only one layer for example introduces this mathematical model. Under the continuous influence of a source with an extracted rate qm-f, the panel solution at position (x, y, z) can be obtained as presented in Supplementary Material 2, and expressed by equation (3)：where pi is initial formation pressure, P(x, y, z) is the pressure at the affected location, B is the oil formation volume factor, ϕm is reservoir porosity, ctm is reservoir total compressibility, α is diffusivity defined as α = km/μϕmctm, x0, y0, and z0 are the position coordinates of a continuous source, and ze is the thickness of a certain layer. Equation (3) only considers the influence of a single fracture element, and then we can calculate the pressure responses affected by the entire vertical fracture using the Superposition principle. The equation is as follows:where x0i,j, y0i,j and z0i,j are the position coordinates of the continuous source (i, j). For the convenience of calculation, we convert equation (4) into dimensionless form. The rewritten equation is as follows:where the definitions of dimensionless parameters in equation (5) are found in Table 1. Finally, applying equation (5) to Nf fracture elements to build the relationship between Nf dimensionless fracture pressures and Nf dimensionless fluxes from the matrix into the fracture.The dimensionless bottomhole pressure can be solved with the model introduced by Peaceman [17], where the wellbore conductivity is assumed to be infinite. Due to fluid exchange behaviors that can occur in each fracture element in contact with the wellbore, the general equation can be given:where pwD is dimensionless bottomhole pressure, qf-wD is dimensionless flux from the fracture into the wellbore, pnwD is the dimensionless pressure of the wellbore, and reqD is the dimensionless equivalent radius. With the aid of equation (6), we have Nj flow equations to participate in the solution of the system of linear equations. Besides, we take into account the impact of the wellbore storage coefficient in the proposed model [18].where γ and CwD are dimensionless variables, defined in Table 1.In this section, one can find that the 2Nf + Nj + 1 flow equations and the 2Nf + Nj + 1 unknowns (Nf dimensionless fracture pressure, Nf dimensionless flux from the matrix to the fracture, Nj dimensionless flux from the fracture to the wellbore and one dimensionless bottomhole pressure) have been given. As such, we can readily solve the system of linear equations using the Gaussian elimination method.To prove the correctness of the proposed model considering finite fracture conductivity and flux variations along the fracture, we compare the calculation results from the proposed model, the analytical model by Bennett et al. [19], and the commercial software (Eclipse). It is worth noting that the flow in the reservoir parallel to the fracture face is ignored in the model of Bennett et al. [19], which means the flow behavior occurring in the reservoir is two-dimensional. Figure 3 presents a top view and a side view of the physical model built in Eclipse, where the red grid represents fracture by setting higher permeability. Significantly, the local grid system belonging to the fracture area is refined to better match the actual results. The reservoir model is discretized into 201 × 201 × 21 grids, each 10 × 10 × 1 m in size. A local grid system of 2 × 1 × 21 grids is further refined into 20 × 5 × 21 grids, among which the width of the fracture is 0.2 m. During the method validation, the following parameters of the matrix, fracture, and fluid are used: layer 1 (ctm1 = 1.2 × 10−5 MPa−1, ϕm1 = 0.1, h1 = 10m, and km1 = 0.01 md), layer 2 (ctm2 = 1.5×10−5 MPa−1, ϕm2 = 0.3, h2 = 10 m, and km2 = 0.03 md), kf = 1 × 104 md, Cw = 0 bbl/psi, Cf = 200 md∙m, ctf = 1.3 × 10−5 MPa−1, ϕf = 0.2, pi = 30 MPa, qw = 0.1 m3/d, μ = 1 mPa∙s, B = 0.985, rw = 0.05 m, xf = 20 m, and zf = 20 m. As shown in Figure 4, the pressure drops and pressure derivatives of the proposed model coincide well with those calculated by commercial software, which proves the validity of the proposed model. However, the calculation results from the model of Bennett et al. [19] are markedly different from the other two, which demonstrates the necessity to consider nonuniform flux distribution.Using the proposed method, a two-layer reservoir is easily constructed to identify the flow regimes observed during the production of an internal well. Subsequently, sensitivity analyses of several parameters on the production dynamic characteristics are performed, including the wellbore storage factor, storability ratio, permeability ratio, fracture height, fracture conductivity, and permeability anisotropy. The detailed parameters used in this section are as follows: xeD = 2, zeD = 2, wD = 0.02, xfD = 2, HD = 2, CfD = 1 × 104, Cs = 0.9630, CwD = 1 × 104, B = 0.985, ω1 = 0.5, ω2 = 1.5, and λ = 1.The curves of dimensionless pressure drop and pressure derivative are plotted in log-log coordinates as shown in Figure 5, in which we can readily recognize four diverse types of flow regimes. (1) Wellbore storage regime [10]. The wellbore storage regime occurs in the early production period, which reflects the magnitude of the fluid elastic energy within the wellbore. The curves of pressure drop and pressure derivative overlap and the value of the slope is 1. (2) Linear flow regime [13]. The linear flow regime is crucial data for determining the fracture conductivity. During this period, a flow behavior that fluid linearly enters the wellbore will happen. The pressure derivative curve has a slope of 1/2. (3) Bilinear flow regime [20]. As the name implies, one more linear flow from the matrix into the fracture happens at the same time as the fracture-wellbore linear flow. This flow regime presents the ability of the fluid within the matrix to flow and is characterized on the pressure derivative curve by a slope of 1/4. (4) Late radial flow regime [9]. During this period, the fluid moves radially into the fracture while the pressure front does not reach the boundary. The characteristic displayed on the pressure derivative curve is a horizontal line (i.e., dPwD/dlntD = 0.5). According to the test data, the boundary properties can be interpreted, including the distance, shape, and orientation.The wellbore storage factor can be used to quantify the magnitude of fluid elastic energy within the wellbore. To investigate the effect of the wellbore storage factors (CwD) on the characteristics of the pressure transient behavior and layered flux contribution, we set CwD = 1 × 103, 1 × 104, and 1 × 105 as shown in Figures 6 and 7. From Figure 6, the calculation results show that the different wellbore storage factors manifest a set of parallel curve clusters to each other. Besides, the wellbore storage coefficient plays a remarkable impact in the early flow stage, and the duration of the wellbore storage regime becomes gradually longer with the increase of CwD value, we can also find that the wellbore storage regime lasts too long to recognize the linear flow regime in the scenario of CwD = 1 × 105. This reflects that the linear flow regime has a high probability of being ignored when the wellbore storage factor is large enough. In Figure 7, one can know that there are significant differences in flux contribution between layers, which is attributed to the layered heterogeneity. The flux contribution of each layer tends to decrease with the increase of the wellbore storage effect. In addition, the flux contributions of the two layers are almost equal early when the CwD with a value of 1 × 105. However, this effect has no effect on the layered flux during the later radial flow regime, which is like the effect on pressure dynamics.According to the definition of the ωn, the lower the storability ratio ω1, the stronger the heterogeneity between the two layers. Figures 8 and 9 depict the pressure dynamics and layered flux curves calculated by the proposed model at different storability ratios. As one can find in Figure 8, the storability ratio has an insignificant effect on the pressure dynamics, and the values of pressure drop and pressure derivative increase slightly as the storability ratio ω1 increases. Therefore, this implies that the multilayer reservoir with stronger heterogeneity will decrease pressure drops and pressure derivatives. As shown in Figure 9, the influence of the storability ratio on flux contribution is mainly concentrated on the linear flow regime and bilinear flow regime. With the increase of the storability ratio (layer 1), the flux contribution of layer 1 gradually increases, while layer 2 is reversed. This reflects that a layer with a larger storability ratio will have a greater flow contribution during mid-production. Besides, note that improving the layer with poor storage capacity only allows it to achieve the same flux contribution as the other layers at the early and middle stages, and the flux contribution will be restored to the proper ratio over time.Permeability is a measure of the ease of passage of fluid through the rock. We investigate the influence on the pressure dynamics and flux contribution curves under three permeability ratios (i.e., χ1 = 0.5, 1.0, and 1.5), as shown in Figures 10 and 11. Compared with Figure 8, we can find that the permeability ratio and storability ratio present the same effect on the pressure drop and pressure derivative curves, that is, the bigger the χ1, the greater the value of pressure drop and pressure derivative. Similarly, this effect on pressure dynamics is not obvious. From Figure 11, one can find that the permeability ratio exhibits a remarkable influence on the layered flux contribution throughout the production period. The flux contribution of layer 1 gradually becomes large as the χ1 increases. In terms of quantitative relations, the flow contribution of different layers is consistent with their permeability ratio. For example, in the scenario of χ1 = 1.5, the permeability of layer 1 becomes three times that of layer 2 and the same is true of the flux contribution ratio between the two layers in the later period. Besides, the flux contribution changes in Figures 9 and 11 before time tD = 1 × 10−2 are very similar, and the difference is that the influence of the permeability ratio persisted until the late radial flow regime. Therefore, this implies that we can really control the flux contribution between layers by changing the permeability ratio.In the production field, the vertical fracture may not fully penetrate the formation, thus it is necessary to investigate the influence of the fracture height on the pressure dynamics and layered flux contribution. Here, we design three cases (zeD = 1.0, zeD = 1.5, and zeD = 2.0) and present the calculation results in Figures 12 and 13. In this study, the position of the fracture is set in the middle of the formation in the vertical direction. In Figure 12, one can observe that the fracture height mainly affects the linear flow regime of the pressure response. A bigger fracture height will result in a smaller value of the pressure drops and pressure derivatives. The reason is that vertical fracture has a larger contact area with the reservoir at a bigger fracture height, resulting in less resistance to flow between the reservoir and the fracture. In Figure 13, we can find that the change of the flux contribution mainly occurs in the early and middle periods. However, the impact of the fracture height is less significant when compared to the above parameters. In addition, with the increase of the fracture height, the flux into the fracture increases at each layer. This is due to that the greater fracture height increases the area of contact with the formation, leading to more channels to allow fluid to enter the fracture.Figures 14 and 15 present the pressure dynamics and flux contribution under different fracture conductivity (i.e., CfD = 1 × 104, 2 × 104, and 4 × 104). In Figure 14, we know that the fracture conductivity exerts an evident impact on the pressure transient behavior at the early and intermediate production. The pressure drops tend to decrease in the face of a higher fracture conductivity, whereas this influence has little effect after the linear flow regime. This is attributed to improved conductivity can reduce the resistance of the fluid flow. In Figure 15, the fracture conductivity affects the flux contribution of each layer in the early and middle periods. The variation characteristic is that the flow contribution of the single layer decreases with the increase of the fracture conductivity, and the effects will disappear at the late radial flow regime.In the pre-assumption of the proposed model, the reservoir is considered to be isotropic. To study the effect of anisotropic permeability on production, we can convert the anisotropic permeability system into an equivalent isotropic permeability system based on the method of Spivey and Lee [21]. In this section, the permeability anisotropy is characterized by a dimensionless parameter Rk (i.e., the ratio of permeability kmx in the x direction to permeability kmy in the y direction). Note that km = (kmxkmy)1/2 is a constant value. Figures 16 and 17 show the pressure data and the layered flux contribution under different permeability ratios, including Rk = 0.1, 1, and 10. It can be seen from Figure 16, that the values of pressure drop and pressure derivative decrease with the permeability ratio increases. This is because the hydraulic fracture propagates along the x direction, and the reservoir permeability along the x direction gradually dominates with the permeability ratio increases. Besides, the results in Figure 17 show that the flux contribution of all layers is monotonically decreased as the permeability ratio increases, whereas monotonically increases over time.In this section, the Mahu-district reservoir is selected as the research area to verify the reliability of the proposed model. Due to the large burial depth and poor permeability, stratified fracturing is widely applied in the Mahu-district reservoir. The proposed model is used to fit the pressure data from a vertically fractured well, which intercepts four oil layers. The well was produced at 17 m3/day, and the other known data of the reservoir and well are summarized in Table 2. Figure 18 shows the measured pressure data and simulated pressure data using the history-matched model. It can be seen that from this figure, the calculated results from the proposed model are in excellent agreement with the field data. Besides, the wellbore after flow and linear flow can be clearly observed in Figure 18. With the help of the history-matching work, the unknown values of the parameters can be determined (as shown in Table 3).A novel semi-analytical method is established to investigate the production dynamic characteristics of a fractured well in multilayer reservoirs. The accuracy of the proposed model is verified with commercial software (Eclipse). Using this model, we successfully differentiate four flow regimes during production. In addition, sensitivity analyses are performed to explore the influence of several parameters on pressure dynamics and flux contributions. The conclusions of this study are listed as follows.The flow regimes that can be distinguished under given reservoir properties and production conditions are as follows: wellbore storage regime, linear flow regime, bilinear flow regime, and late radial flow regime.The presence of a larger wellbore storage factor results in an extended duration of the wellbore storage regime. In addition, the pressure derivative curve is considerably influenced by the value of CwD = 1 × 105, thereby masking the linear flow regime.The storability ratio and permeability ratio exhibit significant influence on the flux contribution of layers. However, it is observed that the permeability ratio can affect the flux contribution of the later period, and the permeability ratio between layers is consistent with the flux contribution ratio between layers.The pressure drop tends to decrease with an increase in fracture height, which indicates a larger fracture can reduce flow resistance. The flux contribution of layers exhibits an increase as the fracture height increases. The reason is that the larger contact area with the matrix.A larger fracture conductivity, as well as a larger permeability ratio, leads to a lower pressure drop.The flow regimes that can be distinguished under given reservoir properties and production conditions are as follows: wellbore storage regime, linear flow regime, bilinear flow regime, and late radial flow regime.The presence of a larger wellbore storage factor results in an extended duration of the wellbore storage regime. In addition, the pressure derivative curve is considerably influenced by the value of CwD = 1 × 105, thereby masking the linear flow regime.The storability ratio and permeability ratio exhibit significant influence on the flux contribution of layers. However, it is observed that the permeability ratio can affect the flux contribution of the later period, and the permeability ratio between layers is consistent with the flux contribution ratio between layers.The pressure drop tends to decrease with an increase in fracture height, which indicates a larger fracture can reduce flow resistance. The flux contribution of layers exhibits an increase as the fracture height increases. The reason is that the larger contact area with the matrix.A larger fracture conductivity, as well as a larger permeability ratio, leads to a lower pressure drop.The data that support the findings of this study are available from the corresponding author.The authors declare that they have no known commercial or associative interest that might influence the work reported in this paper.The authors would like to thank the National Natural Science Foundation of China (No.52104043) for the financial support.Supplementary file S1: Numerical formations of the oil flow in the fracture system.Supplementary file S2: Analytical formations of the oil flow in the matrix system.