Productivity Prediction Model for Stimulated Reservoir Volume Fracturing in Tight Glutenite Reservoir Considering Fluid-Solid Coupling

At present, the main development mode of “horizontal well + volume fracturing” is adopted in a tight glutenite reservoir. Due to the existence of conglomerate, the seepage characteristics are more complex, and the production capacity after volume fracturing is difficult to predict. In order to solve this problem, a dual media unstable seepage model was established for matrix seepage and discrete fracture network seepage when considering the trigger pressure gradient. By considering the Poisson's ratio of stress sensitivity in an innovative way, the coupling model of permeability and stress is improved, and the production prediction model of volume fracturing horizontal well in a tight glutenite reservoir based on the fluid-solid coupling effect is formed. The finite element method is used to numerically solve the model, and the fitting verification of the model is carried out; the impact of stress sensitivity, producing pressure difference, fracture length, number of fracture clusters and fracture flow capacity conductivity on productivity is analyzed, which has certain guiding significance for the efficient development of volume fracturing in a tight glutenite reservoir.


INTRODUCTION
Horizontal well productivity prediction is very important in tight reservoir development. Especially for reservoirs with ultra-low permeability and tight reservoirs, unstable well-production systems may increase the difficulty of well production prediction (Ren et al., 2019;Zhang et al., 2019c). Friehauf et al. (2009) developed a new productivity prediction model that calculates the productivity of a hydraulically fractured well, taking into account the effect of damage of the fracture face from fluid leak-off. Results of the new productivity prediction model are compared with McGuire, Prats and Raymond models. The existing models assume either elliptical or radial flow around the well with permeability varying azimuthally. He found that there were significant differences in the calculated well productivity, which indicate that earlier assumptions made about the flow geometry lead to a serious overestimation of the well PI. The new productivity prediction model can be used to quickly calculate the productivity of wells that have both a finite-conductivity fracture and damage in the invaded zone. Considering the characteristics of unconventional oil and gas seepage, the trilinear flow model was first adopted in the study of volume fracturing horizontal well seepage. The pressure and production dynamics could be analyzed based on the trilinear analytical model, which was adopted by Brown et al. (2011). The two-hole three-region composite analytical model of volume fracturing horizontal wells was established by Brohi et al. (2011). The model described an oil reservoir in the inner region with dual porosity media, oil reservoir in the exterior region with single porosity media. The model was used for the study on tight gas reservoirs, shale gas reservoirs and tight oil reservoirs. Sennhauser et al. (2011) established the multistage fracturing horizontal well productivity model. The model can accurately describe the reservoir in terms of geology, fluid characteristics and pressure profile. It is also demonstrated that the overall productivity of multistage fracturing horizontal wells is determined by a number of parameters, including well spacing, fracture spacing, fracture size and characteristics, and the location of the first and last fractures. The model was applied to the study of tight gas reservoirs, shale gas reservoirs and tight oil reservoirs. Based on the seepage characteristics of porous media in unconventional reservoirs and geological data. Na Zhang and Abushaikha (2019) presented a new fully-implicit mimetic finite difference method (FDM) to simulate reservoir fluid flow. It had been applied to many fields due to its local conservativeness and applicability of any shape of polygon. The principle of the MFD method and the corresponding numerical formula for discrete fracture model are described in detail. The new fully-implicit mimetic finite difference method is tested through some examples to show the accuracy and robustness and the model can be used to predict the well production. Jackson et al. (2013) presented a new approach to simulate fluid flow by adopting Control-Volume-Finite-Element Method (CVFM). The new approach disposes of the pillar-grid concept that has persisted since reservoir simulation began. The new approach promotes the representation of multi-scale geological heterogeneity and the prediction of flow through that heterogeneity significantly. Multiphase flow is simulated using a novel mixed finite element formulation centered on a new family of tetrahedral element types, PN(DG) − PN + 1, which has a discontinuous Nthorder polynomial representation for velocity and a continuous (order N + 1) representation for pressure. The new approach preserves key flow features associated with realistic geological features that are usually lost. Li and Li. (2012) derived a new production prediction model theoretically with the pressure sensitivity of permeability being considered, and used the production data from a low permeability oil field to test and verify the model. The pressure sensitivity coefficient of permeability has been calculated by using the new model with the field data. He founded that the permeability near the well bottom decreased significantly because of the drop in pressure in low permeability reservoirs. An obvious permeability decline funnel could be formed even if the formation was homogeneous before development It was found that the productivity index is no longer a constant in low permeability reservoirs with serious pressure sensitivity of permeability. The pressure sensitivity of permeability should be considered when low permeability reservoirs are being developed. Otherwise, the production will be greatly overestimated. Doe et al. (2013) carried out discretization of complex fracture network in the fractured region, established the productivity model of multiple porosity media, and calculated the contribution rate of natural fractures in the productivity. Cao et al. (2015) established a full analytical mathematical productivity model considering friction, and compared and analyzed the difference between the analytical method and the finite-difference solution results, proving that the model is more accurate in spatial discretization and other aspects . Taking the tight reservoir of Daqing Oilfield as the study object, Hu et al. (2017) established a numerical model of production capacity based on the theory of unsteady seepage mechanics and superposition principle, which was used to calculate and analyze the production capacity of fracturing horizontal wells. Yang et al. (2019) and Zhang et al. (2019a) established the productivity model of the shale reservoir to optimize the fracture spacing, considering the linear flow in the reservoir and the linear flow coupling model in hydraulic fractures, and proved that the productivity of multifracture horizontal wells are inversely proportional to the fracture spacing. The maximum capacity after volume fracturing can be achieved by reducing the fracture spacing. The production capacity of the production well can be doubled by reducing the cluster spacing from 70 to 15 ft.
Compared with the Poisson's ratio of conventional tight sandstone, the Poisson's ratio of tight glutenite reservoir, a mechanic parameter of stone, is greatly different. In this paper, a productivity prediction model for volume fracturing horizontal well considering fluid-solid coupling was established, combining the principle of effective stress and constitutive relation of rock skeleton, while fully taking into account the fluid-solid coupling effect of the tight reservoir; the impact of Poisson's ratio of stress sensitivity and lithological characteristic parameters on the fluid-solid coupling effect was considered in an innovative manner, and a fluid-solid full coupling productivity model suitable for the characteristics of glutenite oil reservoir was established. The accuracy and reliability of the model was validated by comparing the mathematical model of stress field and seepage field of volume fracturing horizontal well in a tight reservoir with the actual production well (Zhang et al., 2019b).

Geometric Model of Volume Fracturing Horizontal Well
In the "Five-region Model" proposed by Stalgorova and Mattar. (2012), and the "Composite Flow Model" proposed by Su Yuliang et al., it is believed that the main fracture and natural fracture are formed by fracturing form complex fracture clusters after volume fracturing of horizontal wells, and there are non-fractured regions between the clusters. The reservoir area of volume fracturing horizontal wells are composed of three parts: reservoir matrix, natural fracture (unmodified) region, artificial fracture and natural fracture interlaced fracture network region. Based on the above conclusions, a volume fracturing horizontal well model for a tight glutenite reservoir is established, as shown in Figure 1.
The formation and fluid meet the following conditions: 1) thickness of the reservoir is h, with outer boundary closed, and with natural fracture; 2) rocks and fluids are slightly compressible; 3) the impact of gravity and temperature changes are not considered for the seepage process.

Mathematical Model for Seepage of Volume Fracturing Horizontal Well
Mathematical Model for Seepage of Matrix-Natural Fracture System The fluid motion equation in a matrix system considering the trigger pressure gradient: Equation for state of rock skeleton and fluid in matrix system Equation for continuity of a single phase compressible fluid Eqs (1) and (2) were substituted into the continuity Eq. (3): Since the ∇P m is very small, Eq.
(3) can be reduced to the following equation.
Eq. (6) is the equation for the seepage governing differential of single-phase micro-compressible fluid in bedrock when the trigger pressure gradient is considered. Similarly, the governing differential equation of the natural fracture system can be obtained as follows.
The initial and boundary conditions:

Seepage Model of Fracture Network System
The fluid motion equation: The equation for state: The equation for continuity Similarly, the seepage governing differential equation: The initial and boundary conditions P f x, y, z, t 0 P i P f x, y, z, t P m x, y, z, t P n x, y, z, t

Definite Conditions for Solution of Seepage Field
In order to solve the nonlinear equations, some fixed value conditions are set. These conditions are definite conditions. For example, the fixed pressure condition is that the well is produced under a certain flow pressure. For the mathematical model for the seepage field of volume fracturing horizontal wells in tight reservoirs, the definite conditions mainly include the initial formation pressure conditions, the definite pressure production. Initial conditions: Boundary conditions:

Model for Stress Field of Volume Fracturing Horizontal Well in Tight Reservoir
Assumptions: 1) The rock is the deformation of porous media; 2) The rock particles are incompressible, while particle pores can be compressed; 3) The impact of temperature change on rock deformation is not considered; 4) Rock deformation is a small elastoplastic deformation; 5) The compression coefficient of pore changes (Jianjun et al., 2002;Liu et al., 2018;Zhao and Du, 2019).

Equation for Constitutive Relation
Geometric Equation

Differential Equation for Stress Equilibrium
Before the reservoir is developed, any unit (infinitely small) in the reservoir is in static equilibrium. When the reservoir is developed, the pore pressure of the reservoir will change. Meanwhile the effective stress also changes. This change can be expressed by the stress equilibrium differential equation. A hexahedral differential element is used to represent the stress state in the space, as shown in Figure 2 and the length of each side is dx, dy, dz respectively. The equilibrium equation of the element in the x direction: So if we simplify Eq. (18), we get Eq. (19).
Similarly, the differential equation for equilibrium in the y and z directions:

Definite Conditions for Solution of Stress Field
For the mathematical model for stress field of horizontal well in a tight reservoir, its definite conditions can be divided into initial conditions and boundary conditions, and the boundary conditions can be divided into stress boundary conditions and displacement boundary conditions (Zhang et al., 2018). The initial conditions for the stress field are consistent with those for the model for seepage field. Stress field displacement boundary conditions: Stress field stress boundary conditions: Substitute the data into Eqs 24-27: Coupling Relation of Porosity Ran and Li. (1997) derived the porosity strain model through the definition of porosity. The volumetric strain parameters themselves imply the mechanical properties of rocks, and the derived porosity coupling model is as follows:

Finite Element Solution of Seepage Model
The whole reservoir is considered as a continuous calculation domain, including porous media seepage considering the trigger pressure gradient and the fracture network seepage of the volume fracturing region, and fracture fluid flow within the volume fracture system (Zhang et al., 2020). For fracture fluid seepage in the fracture network system, the actual a three-dimensional fracture is equivalent to two-dimensional fracture with certain fracture openness with a discrete fracture network model. So, the whole calculation domain governing equation Ω i,m of threedimensional integral form can be reduced as.
The approximate expression of the pressure field function is obtained by taking any element in the system: Eq. (30) is substituted into Eq. (6) to obtain the expression of characteristic matrix of element Ω i,m ∭ Ωi,mn Similarly, the element characteristic matrix expression of natural fracture seepage model Ω i,n : The element characteristic matrix expression of network fracture seepage Ω i,f in volume fracturing region: Assuming that the initial reservoir flow occurs in the natural fracture system Ω i,n and flows into the fracture network system Ω i,f , the solution can be obtained by the implicit backward difference scheme, and finite element closed equations with unknown quantities in the next time step Ω i,f can be obtained: The pressure at the next time step of the matrix system is calculated according to the matrix natural fracture equilibrium equation of the entire reservoir: The whole continuous calculation domain can be solved through Eqs (34) and (35).

Finite Element Solution for Fluid-Solid Coupling Productivity Model
Differential equation for stress field: Boundary conditions: The unknown function U is the field function that needs to be solved. A and B in Eqs (36) and (37) are differential operators for independent variables such as spatial coordinates and time. The equivalent integral form of differential Eq. (36) and boundary conditions Eq. (37) can be written by the Galerkin equation as follows: According to the principle of effective stress, the total stress of reservoir rock is composed of effective stress and pore fluid pressure: It is assumed that the displacement vector at any point in the calculation domain satisfies the following equation: The above equation can be expressed as: The following equation can be obtained by the combination of Eqs (17) and (41): By substituting Eqs (19)- (23) and Eq. (36) into Eq. (38), the following equation can be obtained: The above equation is integrated by parts and substituted into Eq. (42), Eq. (17) and the elastic constitutive matrix of rock to obtain the overall equilibrium equation of the fluid-solid fully coupled stress field: The Galerkin weighted residual method is applied to the seepage field model and the following equations can be obtained: To sum up, the matrix form of the equilibrium equations for the fluid-solid coupling mathematical model of porous media based on the full-coupling finite element method can be obtained as follows: Eq. (51) is substituted into COMSOL and the result is obtained.

Comparison of Actual and Model Production
After establishing and solving the fluid-solid coupling model, the required parameters of the model are input into COMSOL ( Table 1). Then the trend of the model simulated oil rate is obtained by using the full coupling of finite elements. Combined with the actual oil rate changes of production well, the trend of model simulated oil rate can be analyzed as shown in Figure 3. It can be seen from Figure 3 that the fitting error between the actual daily production and the simulation result of fluid-solid coupling model is small.  Frontiers in Energy Research | www.frontiersin.org November 2020 | Volume 8 | Article 573817 7 distribution at 0, 300 and 600 days of production, respectively. As can be seen from the Figures, formation pressure gradually decreases in the horizontal well simulation production process, and the pressure change range is as follows: area near the wellbore of horizontal well ≥ volume fracturing region ≥ unmodified region of matrix. Globally, the reservoir fluid flow first occurs around the wellbore with a fast flow speed. After the fluid flows away, the supply fluid of the volume fracturing region and matrix flows into the wellbore from the fracturing region and matrix flow into the wellbore. The fluid of the entire formation flows around the wellbore of the horizontal well elliptically.

Sensitivity Analysis
By sensitivity analysis of different parameters, the impact of each factor on production can be understood. These parameters mainly include stress sensitivity, producing pressure, fracture half length, number of fracture clusters and fracture flow capacity.

Stress Sensitivity
The fluid-solid coupling effect is mainly due to the change of reservoir pressure in the production process, which leads to the change of pore and permeability. And the change of permeability is obvious. Based on the rock mechanics characteristics of the glutenite, the stress sensitivity coefficient is a parameter to reflect the influence of stress change on permeability. The flow pressure is set to 26 MPa and other parameters remained unchanged by single factor analysis. The variation of the oil rate under different stress sensitivity coefficients is simulated.
As can be seen from Figures 7, 8, with the increase of the stress sensitivity coefficient, the production of the horizontal well is  reduced. Under the effect of stress, the initial permeability of the reservoir, especially the volume fracturing region, is high, and the permeability loss with the change of stress is large, resulting in the decrease of production. In the actual production process, the control of the producing pressure should be considered to reduce the fluid-solid coupling effect on the productivity.

Producing Pressure Difference
The change of pressure in the production process will also lead to the change of pore and permeability of reservoir. The fluid-solid coupling can be reflected by controlling the producing pressure difference of the model. The stress sensitivity coefficient is set to 0.06 MPa −1 and other parameters remained unchanged by single factor analysis. The variation of oil rate under varying producing pressure differences is simulated.
As can be seen from Figures 9, 10, the daily production of volume fracturing horizontal well increases with the increase of producing pressure difference. The increase of producing pressure difference leads to a large drop in formation pressure, which will obviously affect the loss of formation porosity and  permeability. Increasing the producing pressure difference will increase the cumulative production of the horizontal well, but the increasing rate will decrease. After increasing the producing pressure difference, the reservoir's available range increases, leading to higher production. However, it also enhances the fluid-solid coupling effect, resulting in a decrease in fluid flow capacity, which has some impact on the actual production of the horizontal well.

Fracture Half Length
Fracture half length is one of the important parameters to control the area of SRV zone. By extending the fracture half length, the SRV zone of the reservoir can be expanded effectively. Therefore, fracture half-length has a great influence on horizontal well production. The variation of the oil rate under different fracture half lengths is simulated by single factor analysis.  As can be seen from Figures 11, 12, the increase of fracture half length has a small impact on the early productivity of the horizontal well. With the increase of fracture half-length, the growth trend of cumulative production is obvious. However, although the fracture half-length increases in a linear manner, the growth rate of cumulative production decreases to some extent. At the same time, considering the economic cost of on-site fracturing, there is some room for optimization of fracture half-length.

Number of Fracture Clusters
The number of fracture clusters could control the area of SRV zone. Therefore, the number of fracture clusters has a great influence on horizontal well production. The variation of oil the rate under different numbers of fracture clusters is simulated by single factor analysis.
As can be seen from Figures 13, 14, with the increase of fracture clusters, the growth trend of cumulative production is obvious. The increasing of fracture clusters has a great impact on the production of  the horizontal well in the early stage and a small impact in the later stage. In the early production of the horizontal well, the oil supply is mainly from the volume fracturing region to the horizontal wellbore, and the oil supply is from the reservoir matrix region in the later stage. Due to the unmodified reservoir matrix, the seepage condition is poor and the oil supply capacity is relatively weak, resulting in little difference in the daily production at the later stage.

Fracture Flow Capacity
Fracture flow capacity is another important parameter to control horizontal well production. The variation of the oil rate under different fracture flow capacities is simulated by single factor analysis.
As can be seen from Figures 15, 16, the fracture flow capacity has a great impact on the early production stage of the horizontal well and a small impact on the later stage. In the early stage of horizontal well production, the wellbore of the horizontal well is mainly supplied from the volume fracturing region, which results in a high fracture flow capacity and fast fluid flow in the effective modification area. In the later stage, oil is supplied from the reservoir matrix region, which is not modified, has low permeability and relatively weak oil supply capacity, resulting  in little difference in daily production in the later stage. However, although the fracture flow capacity increased in a linear manner, the growth rate in cumulative production decreased.

CONCLUSION
According to the seepage characteristics, the mathematical model of multiple media seepage in volume fracturing horizontal wells is established. The model takes into account the dual media composed of the fracture grid and the matrix formed by the coupling of artificial fracture and natural fracture, and a seepage model for the matrix region considering the trigger pressure gradient is established. In addition, the fluid-solid coupling effect in the tight reservoir is considered in the model, and the basic mathematical model of the stress field is established based on the principle of effective stress and the constitutive relation of rock skeleton. In addition, the mathematical model of stress field-seepage field for volume fracturing horizontal wells in tight reservoir is formed by combining with the seepage model. Finally, based on the rock mechanics characteristics of glutenite and taking into account the impact of Young's modulus, Poisson's ratio, initial permeability value and other parameters, the permeability stress crossing dynamic model is improved, and the fluidsolid coupling productivity model for volume fracturing horizontal wells is formed.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
All the authors have contributed to the conception and design of the productivity model. LT was mainly responsible for the construction of the capacity model, XC and PW were responsible for the solution and validation of the capacity model, and HW was responsible for the sensitivity analysis of the capacity model. The first draft was written by XC, and was revised and confirmed by all the authors. GLOSSARY K m matrix permeability, 10 −3 μm 2 v m fluid flow rate in matrix, 10 −3 m/s μ fluid viscosity, mPa · s ∇P m the pressure gradient of the pore system in bedrock, MPa/m λ trigger pressure gradient, MPa/m ϕ matrix porosity, % ϕ 0 initial formation porosity, % C L fluid compressibility, MPa −1 C p pore compressibility, MPa −1 P i initial reservoir pressure, MPa ρ fluid density, kg/m 3 ρ 0 initial fluid density, kg/m 3 Q m cross flow rate, kg/m 3 ·s α = L 2 m /12 form factor L m matrix rock size, m C m matrix composite compressibility, MPa −1 P n natural fracture system pressure, MPa P m matrix system pressure, MPa K n natural fracture permeability, 10 −3 μm 2 C n natural fracture compressibility, MPa −1 ϕ n natural fracture porosity K f artificial fracture permeability, 10 −3 μm 2 C f artificial fracture compressibility, MPa −1 ϕ f artificial fracture porosity p w bottom hole flowing pressure, MPa dσ ij effective stress increment D ijkl elastic-plastic coefficient matrix tensor dε kl strain increment ε ij strain tensor μ displacement component f i volume force on one side σ normal stress τ shearing strength τ xy the component of stress in the y direction on the x plane α effective stress sensitivity coefficient, MPa −1 E Young's modulus, GPa υ Poisson's ratio α E relationship between Young's modulus and effective stress sensitivity coefficient α υ relationship between Poisson's ratio and effective stress sensitivity coefficient α K relationship between initial permeability value and effective stress sensitivity coefficient c 1 Young's modulus influence coefficient c 2 Poisson's ratio influence coefficient c 3 initial permeability value influence coefficient ε v volumetric strain N i,mn system element basis function matrix of Ω m,n ∇N T i,mn , ∇N i,mn transpose of matrix P i,m element node pressure matrix of Ω m V any function on the boundary σ total stress of reservoir rock σ 9 effective stress M System mass matrix F Γ load vector C System damping matrix L horizontal well length, m σ max maximum principal stress, MPa σ min minor principal stress, MPa Frontiers in Energy Research | www.frontiersin.org November 2020 | Volume 8 | Article 573817