Nonlinear Seepage Mathematical Model of Fractured Tight Stress Sensitive Reservoir and Its Application

A tight reservoir can only have industrial production value after fracturing. After volume fracturing, there are artificial fracturing zones and fracturing unmodified zones composed of a matrix, microfractures, and main fractures. There are also pressure-sensitive effects and nonlinear seepage in the matrix. There is an unsteady channeling flow between the matrix and microfractures and main fractures and a threshold pressure gradient in the matrix. The fluid flow between the three media and the transfer flow between the media are very complex. Aiming at the problem of poor simulation accuracy because the above characteristics of tight reservoir fracturing development are not considered in the current commercial numerical simulation software Eclipse, CMG, and VIP, the triple medium coupling nonlinear mathematical model of tight reservoir fracturing development matrix/microfractures/main fractures was established in this study, and the discrete mathematical model was constructed using the full implicit numerical solution method. The model considered the effects of threshold fracturing gradient and nonlinear seepage on the reverse side of seepage. In terms of pressure-sensitive effect, the latest generalized pressure-sensitive effect calculation model, considering rock physical parameters established by authors, was adopted. In terms of matrix/fracture transfer flow, the nonlinear matrix/fracture transfer flow model considering threshold pressure gradient and pressure-sensitive effect established by authors was adopted. Finally, a tight reservoir numerical simulation software was developed based on a triple medium coupling nonlinear model. The software was used to simulate the development effect of turbidite tight reservoir in Wang 587 Block of Jiyang depression, Shengli Oilfield. The results showed that the fitting accuracy of the new model was 13.3% higher than that of Eclipse, and the calculation results were more in line with the production practice of tight reservoirs. The new model established an important theoretical basis for scientifically and effectively guiding the development of tight reservoirs.


INTRODUCTION
Reservoir numerical simulation technology is an effective means to guide development scheme design, development adjustment, and scheme prediction and provide a theoretical basis for oilfield development scheme preparation and policy-making (Chengzao et al., 2012a;Caineng et al., 2012). A tight reservoir is usually defined as a reservoir with gas-measured permeability less than 1.0 × 10 -3 μm 2 . It has also been suggested that this limit is not very strict because a tight reservoir is defined as a reservoir with gasmeasured permeability less than 2.0 × 10 -3 μm 2 based on the description of the National Energy Administration of China (Caineng et al., 2013;Fuxi et al., 2014). Taking the Jiyang Depression of Shengli Oilfield studied in this article as an example, the field classifies glutenite, turbidite, and beach bar sand reservoirs with gas-measured permeability less than 3.0 × 10 −3 μm 2 as tight reservoirs for development (Chengzao et al., 2012b;Zhenglian et al., 2012). Due to the differences in reservoir genesis of the three types of tight reservoirs, different types of tight reservoirs are developed in different ways. At present, effective development technology of three-dimensional injection and production of sand gravel "multi-storey" multistage fracturing horizontal wells, effective development technology of segmented fracturing fracture network adaptation of turbidite vertical wells, and long fracture fracturing development technology of beach bar sand vertical wells have been developed. However, due to the lack of special numerical simulation software for tight reservoirs, the development plan lacks a theoretical basis, and individual production blocks have problems such as large well spacing, low formation pressure level, and slow oil production rate. Special tight reservoir numerical simulation software is required for development simulation (Clarkson et al., 2013;Zhongxiang et al., 2016). In order to introduce the development characteristics of tight reservoir artificial fracturing zones, this study investigated the research status of pressure-sensitive effect, nonlinear seepage, nonlinear transfer flow, and tight reservoir numerical simulation.

Pressure-Sensitive Effect
In the study of pressure-sensitive effects in tight reservoirs, scholars currently mainly use experimental data fitting methods when studying pressure-sensitive effects. For experimental research methods, Holt et al. and Tian et al. experimentally studied the stress sensitivity of matrix cores and fractured natural cores and gave an exponential empirical model (Holt, 2005). Nie et al. investigated the effect of time effect on the stress sensitivity of tight gas reservoirs and established an empirical model of stress sensitivity considering the time factor (Xiangrong et al., 2018). Yang et al. combined CT structural scanning experiments to quantitatively evaluate the variation of pore and throat radii with surrounding pressure (Yang et al., 2017). Gu et al. experimentally analyzed the main controlling factors of stress sensitivity in tight sandstone in terms of material composition, fracture development, and pore structure. However, they did not give a quantitative model (Yang et al., 2018). Liu et al., based on experimental results, considered the relationship between effective stress sensitivity coefficient and Young's modulus and established an exponential stress sensitivity model (Liu et al., 2018). Yang et al. analyzed the variation of pore structure parameters with effective stress from the microscopic scale by combining CT structural scanning techniques and gave empirical equations for the stress sensitivity of porosity and permeability (Yongfei et al., 2016). Tian et al. quantified the relationship between strain and effective stress based on experimental results and considering the effects of boundary layer and critical throat radii (Tian et al., 2015). Lei et al. applied an exponential-type function to quantify the stress sensitivity of shale reservoirs considering the effect of the initial permeability value (Hao et al., 2019). Shangxian Yin and Jones et al. proposed empirical models of permeability and effective stress in different functional forms. However, various empirical regression models were only applicable to a single sedimentary rock type and cannot be extended to different sedimentary rock reservoirs (Jones, 1975;David et al., 1994;Dong et al., 2010).
In the study of the pressure-sensitive effect of microfractures in tight reservoirs, Yanlong Zhao, Zijun Zheng, and others carried out a fracture permeability stress sensitivity experiment and analyzed the fracture structure characteristics by scanning electron microscope technology. They proposed that fracture permeability and effective stress show a negative exponential relationship (Zijun and Cheng, 2018;Yanlong et al., 2019). Lanqing Fu and others studied the stress sensitivity characteristics of tight cores with fractures under the experimental conditions of constant displacement pressure difference and constant confining pressure, respectively. They proposed that there was an inflection point in the fracture stress sensitivity curve, which was before the effective stress reached the inflection point, the permeability decreased sharply, the effective stress exceeded the inflection point, and the decline range of permeability tended to be gentle (Lanqing, 2016). Pingliang Fang studied the permeability loss of micro and nano fractures under variable pore pressure and variable confining pressure. The permeability loss rate exceeded 70%, and the damage degree of stress sensitivity was "strong" (Jianwei et al., 2019). Jianwei Feng and others proposed that, with the increase in effective stress, the fracture permeability gradually decreased and there was an "inflection point." The reason was that brittle fracture occurred after the volume expansion of brittle minerals underground conditions, the connecting dominant channel was damaged, and the fracture system was complex. Then, the mechanical relationship model between crack opening and effective stress was deduced by the mathematical method. It was clear that the fracture opening conformed to the nonlinear variation law with confining pressure and effective stress . Wei Tian et al. studied the stress sensitivity of the reticulated fracture core using a fully digital servo controller, obtained the conclusion that the shape of the stress sensitivity curve presents an "L" shape, and divided the curve into two stages. The first stage was pseudoplastic deformation, and the permeability decreased greatly. The second stage was dominated by elastic deformation. As the rock became an isobaric body, the decline range of permeability became smaller (Pingliang et al., 2017a). Fraser Harris et al. studied the effect of fracture radial displacement and shear displacement on microfracture permeability by changing the size and direction of the radial stress field (Fraser Harris et al., 2020). Liang et al. established exponential microfracture pressure-sensitive model and cubic polynomial variation law of main fracture conductivity through microfracture and main fracture stress sensitivity experiments. They proposed that the influence of the microfracture pressure-sensitive effect on fluid flow was greater than that of the main fracture conductivity loss (Bin et al., 2016). He et al. studied the pressure-sensitive effect of microfractures through experiments and divided its stress sensitivity into medium to strong according to permeability loss. Combined with CT scanning imaging experiments, they studied the influence of the fracture surface on the pressuresensitive effect of microfractures (Jingang et al., 2013). Haiyong Zhang et al. prepared the thermal stress microfracture dense core by heating in a microwave oven to approximately simulate the microfracture formed in the actual diagenesis and sedimentation process. Through the indoor experiment of stress sensitivity of the microfracture core, the stress sensitivity was evaluated as medium to strong, but the pressure-sensitive model was not given (Hai Yong et al., 2015).

Nonlinear Seepage
In the study of nonlinear seepage equations in tight reservoirs, Yanzhang Huang proposed a boundary layer theory that when the fluid flowed in micro-and nano-pore throats, there existed both bulk phase layers and boundary layers. When the driving pressure difference was low, the fluid in the boundary layer did not flow. With the gradual increase in pressure difference, the thickness of the boundary layer gradually decreased, which was characterized by nonlinear seepage under small pressure difference and linear seepage characteristics under large differential pressure (Shiming, 2019), as shown in Figure 1. In addition, scholars have verified this theory by conducting physical simulation experiments. Through core displacement experiments, Changjun Zhu et al. verified that the smaller the permeability, the greater the proportion of nano-throats in the total pore volume because with the decrease of permeability, the starting pressure gradient increased, the surface tension increased, and the nonlinear seepage became more obvious (Changjun, 2002). Pu Zhang et al. applied experimental methods to study the boundary layer on the seepage pattern and concluded that the smaller the radii of the throats, the greater the influence of the boundary layer on nonlinear seepage (Pu et al., 2008). However, in practical applications, it was difficult to accurately measure the boundary layer thickness by experimental means. Therefore, the seepage equations based on the boundary layer theory were limited in numerical simulation applications. At present, formula v − K μ (∇P − D P ) with pseudo threshold pressure gradient (D P ) is widely used to describe non-Darcy seepage in tight reservoirs. This formula has two problems: the experimentally determined proposed threshold pressure gradient is significantly high and it cannot describe the nonlinear section of low-velocity seepage.

Nonlinear Transfer Flow
In studying the matrix/fracture transfer flow laws, the most commonly used are Warren and Root and Kazemi formulations. Warren and Root proposed a constant shape factor formulation α 4n(n+2) L 2 . In 1976, Kazemi (Kazemi et al., 1976) applied a finite difference equation to obtain a shape factor expression for each homogeneity in three dimensions α 4( 1 which was the most commonly used dual-medium transfer flow equation for commercial numerical simulation software.
The matrix/fracture transfer flow laws of tight reservoirs differed from those of conventional sandstone reservoirs. The pressure changes of the matrix/fracture system were unsteady and affected by the pressure-sensitive effect, as shown in Figure 2.
The pressure-sensitive effect on the transfer flow was mainly reflected in the following two aspects: first, the pore pressure change and permeability loss caused by the production pressure difference change during the development of tight reservoirs affected the matrix/fracture transfer flow. Second, in the process of transfer flow, the average matrix pressure gradually decreased, and the decrease in pore pressure caused permeability loss, thus affecting the instantaneous transfer flow.  Relatively few studies on non-steady-state transfer flow in tight reservoirs were conducted at home and abroad, and the related studies were mainly adopted theoretical studies, lacking experimental verification. In 2007, Yongming He derived the matrix/fracture transfer flow shape factors for one-dimensional, two-dimensional, and three-dimensional cases, considering the effect of the threshold pressure gradient (Yongming et al., 2017), as shown in Eq. 1. However, when establishing the pressure diffusion equation, the seepage equation was quasilinear, which cannot well describe the low-speed non-Darcy seepage stage of tight reservoirs: ( 1 ) In 2014, Weibing Shu et al. developed a matrix/fracture transfer flow model by fitting experimental data on seepage in extra-low-permeability cores and solved for the transfer flow shape factor (Weibing et al., 2014), as shown in Eq. 2. However, they did not consider the initiating fracture gradients and pressure-sensitive effects: (2) In 2018, Shenglai Yang established the matrix pressure diffusion equation considering the pressure-sensitive effect in tight reservoirs and applied the Boltzmann transformation to solve the transfer flow shape factor (Kazemi and Gilman, 1993), as follows: However, there were still some limitations in its conditions of use. First, a simple exponential pressure-sensitive model was used, which cannot be extended to different types of tight reservoirs. Second, a linear seepage equation was used to establish the matrix pressure diffusion equation, which cannot characterize the non-steady-state pressure diffusion characteristics in the fracturing process and the influence of the threshold pressure gradient.
Through the above research, different scholars have studied the dual medium matrix/fracture transfer flow law and influencing factors from multiple perspectives. However, most of the research results were only applicable to conventional sandstone reservoirs, while not to tight reservoirs with nonlinear seepage characteristics and pressure-sensitive effects.

Numerical Simulation of Tight Reservoirs
In the theoretical research of tight reservoir numerical simulation, Pingliang Fang et al. considered the evolution laws of matrix elastic shrinkage, fracture initiation and closure, and fracture and extension in the fracturing process and production process and established the multi-media dynamic change mechanical mechanism and fluid-solid coupling mathematical model of tight sandstone reservoirs (Pingliang et al., 2017b). Aiming at the problem of repeated fracturing in tight reservoirs, Kang Ma et al. established a two-dimensional two-phase nonlinear seepage mathematical model using the discrete main fracture treatment method (Kang et al., 2015). Dongling Yang established a threedimensional three-phase seepage mathematical model of fractured reservoirs and developed a numerical simulator using the full implicit solution method. Both the seepage equation and matrix/fracture seepage equation are linear (Dongling, 2014). Based on the nonlinear seepage theory, Lu Li established the nonlinear mathematical model of tight reservoirs by considering the end effect of capillary force, steady matrix/fracture transfer flow, and tight reservoir imbibition (Lu, 2020). Considering the start-up pressure gradient, pressure-sensitive effect, and imbibition of tight reservoirs, Hong Li established a linear seepage model and developed a tight reservoir numerical simulator (Hong, 2015). Based on the experimental data, Rongze Yu and Huaqiang Zhang established the nonlinear seepage equation and the threedimensional three-phase nonlinear seepage mathematical model of ultra-low permeability reservoirs (Rongze et al., 2012). Yueli Feng considered the pressure-sensitive equation of fracture directionality, established the mathematical model of oil, gas, and water three-phase seepage, considering the pressuresensitive effect of fractures in different directions, and studied the influence of fracture angle and pressure-sensitive effect on development effect through numerical simulation (Yueli, 2019). Using the analytical solution method and effective numerical calculation method for solving the moving boundary problem in heat transfer, Wenchao Liu conducted accurate numerical simulation research on single-phase or oil, gas, and water three-phase nonlinear seepage in low permeability reservoirs based on the quasilinear motion equation and nonlinear motion equation (Wenchao, 2013). However, the above model had no substantive breakthrough in considering the pressure-sensitive effect, nonlinear seepage, and unsteady matrix/fracture transfer flow.
Based on the above investigation, the commercial numerical simulation software had some defects in simulating tight reservoirs. First, the mathematical model was generally derived based on the Darcy seepage equation, which cannot characterize the nonlinear seepage characteristics of a tight reservoir matrix. Second, the pressure-sensitive model was only a simple exponential function that could not be applied to different types of tight reservoirs. Thirdly, the dual medium steady-state transfer flow model proposed by Warren and Root could not characterize the unsteady transfer flow law of matrix/ fracture in tight reservoirs. Therefore, the matrix/ microfracture/main fracture triple medium mathematical model of tight reservoir fracturing development considering the pressure-sensitive effect, nonlinear seepage, and nonlinear channeling flow was established. Herein, the discrete mathematical model based on the full implicit solution method was constructed, and the special numerical simulation software suitable for tight reservoirs was developed. This laid a foundation for guiding the largescale and efficient development of tight reservoirs.

ESTABLISHMENT OF MATHEMATICAL MODEL
In the field of numerical simulation theory of tight reservoirs, scholars generally believe that there are three media of matrix, microfractures, and main fractures in tight reservoirs, as shown in Figure 3. Considering the influence of the pressure-sensitive effect on the seepage characteristics of different media and the nonlinear characteristics of matrix seepage and matrix/fracture transfer flow, a triple media coupling mathematical model was established in this study.

Characterization of Pressure-Sensitive Effect
In 2020, Kai Liu and Dayin Yin et al. established a theoretical model of the pressure-sensitive effect applicable to different tight reservoirs through theoretical analysis based on experimental data (Kai et al., 2020a), as shown in Eq. 6. The specific idea is to determine the relationship between permeability and pore structure parameters (pore radii and throat radii) based on the Kozeny-Carman equation. Based on the Hertz contact deformation theory, the relationship between pore radii and effective stress was determined. Combined with microscopic experiments, the relationship between throat radii and effective stress was determined. Substituting the relationship between pore radii and throat radii into the permeability expression, the calculation formula of matrix pressure-sensitive effect was obtained: After simplifying the experimental data of rock physical properties of tight reservoirs in Jiyang depression, the formula can be written as ξ m (p m ) could be calculated as the following equation: In 2020, Kai Liu, Daiyin Yin, and others established an empirical model for calculating the pressure-sensitive effect of microfractures based on the pressure-sensitive experiment of microfractures as follows: Nonlinear Seepage Characteristics In 2019, Fuquan Song deduced the nonlinear seepage equation of the tight reservoirs based on the experimental results of nanotube microseepage and the boundary layer theory: After Taylor expansion, the higher-order term was ignored, and the pressure-sensitive calculation model Eq. 6 was introduced. Under the pressure-sensitive condition of the tight reservoirs, the nonlinear seepage equation was Based on the relative permeability curves, the single-phase seepage equation was extended to three-phase flow, and the general equation of seepage for oil, water, and gas can be written as Where p is the oil, water, and gas phases and k rp is the relative permeability of the oil, water, and gas phases.  Figure 4.

Characterization of Nonlinear Transfer Flow
The results of the comparison between the new transfer flow model and the conventional model showed the following: 1) the average error of the 1D and 2D instantaneous transfer flow calculation of the modified transfer flow theory model was 6.6 and 7.8%, which was more than double the accuracy of the classical Warren and Root model and 26.5% higher than that of the Kazemi model. 2) The transfer flow calculation results of the conventional Warren and Root model were always large. The reason for this was that, on the one hand, the transfer flow equation was based on the Darcy seepage, which did not consider the influence of the threshold pressure of the tight reservoirs, and the pressure difference was large. On the other hand, the derivation of this model was mainly used in well-testing theory, and the stage of the well-testing interpretation was usually the stage of pressure recovery instability, so the given empirical equation of the transfer flow equation exaggerates the elastic energy.
3) The transfer flow calculation result of the Kazemi model was low in the early stage and high in the late stage, which was because the transfer flow equation was based on the stable seepage theory, the effects of the tight reservoirs threshold pressure and elastic energy were not considered.
The new transfer flow model was derived based on the nonlinear seepage equation, which considered the threshold pressure and pressure-sensitive effects on the tight reservoirs and was applicable to the study of matrix/fracture transfer flow in the tight reservoirs. In contrast, the conventional Warren and Root model and Kazemi model were derived based on the steady seepage theory and applied to conventional sandstone reservoirs while the application to the tight reservoirs had some limitations.

Mass Conservation Equation
Herein, three-dimensional three-phase mass conservation equations of matrix, microfractures, and main fractures were established by considering the matrix pressure-sensitive effect, microfracture pressure-sensitive effect, matrix nonlinear seepage characteristics, and matrix/fracture nonlinear transfer flow in fracturing development of tight reservoirs.
The matrix mass conservation equation are as follows: Oil phase: Water phase: Gas phase: ϕ m is matrix porosity; γ o , γ w , γ g are the gravity coefficients of the oil, water, and gas phases; D is the depth calculated from a datum, downward is positive, m; S m o , S m w , S m g are the saturation of the oil, water, and gas phases in a matrix; T m o , T m w , T m g are the conductivities of the oil, water, and gas phases flowing from reservoir grid i to adjacent reservoir grid i + 1; and R s is dissolved gas-oil ratio. The conservation equation of microfractures mass is as follows: Oil phase: Water phase: Gas phase: ϕ f is the microfracture porosity; q fF o , q fF w , q fF g are volume transfer flows of the oil, water, and gas phases between microfractures and main fractures, m 3 /s; S f o , S f w , S f g are the saturation of the oil, water, and gas phases in microfractures; T f o , T f w , T f g are the conductivities of the oil, water, and gas phases flowing from microfracture grid i to adjacent grid i + 1 in the microfracture.
The mass conservation equation of main fractures are as follows: Oil phase: Water phase: Gas phase: ϕ F is main fracture porosity; S F o , S F w , S F g are the saturation of the oil, water, and gas phases in the main fractures; q FW o , q FW w , q FW g are volume flows of the water, oil, and gas phases from the main fractures to horizontal wellbores per unit time, m 3 /s; T F o , T F w , T F g are the conductivities of the oil, water, and gas phases flowing from fracture grid i to adjacent fracture grid i + 1 in the main fracture. All fracture grids have flow exchange with adjacent microfracture grids on both sides, and only the main fracture grid passing through the horizontal wellbore has a flow term ofq FW o , q FW w , q FW g , and the transfer flow terms of other fracture grids are zero.

MATHEMATICAL MODEL SOLUTION
There are two kinds of numerical methods for solving multiphase seepage problems: one is the implicit pressure explicit saturation method (IMPES method) and the other is the simultaneous method; that is, by solving simultaneous equations, each independent variable is solved at the same time. When using the IMPES method, the calculation is small and the form is simple. However, for the multi-media numerical simulation of tight reservoirs, the variable permeability coefficient, and pressure are coupled in the seepage equation, and the convergence of the model is poor. The fully implicit simultaneous solution method with high stability and better convergence is adopted in this study.

Fully Implicit Simultaneous Solution Method
The fully implicit simultaneous solution method requires multiple iterations at each time step to solve the multiple coefficient matrix. Through the implicit difference of the mathematical model, the nonlinear difference equations of matrix, microfractures main fractures, and horizontal wellbores were established. After linearization by the fully implicit method, the general expression of the coefficient matrix of each iterative step in each time step was obtained, and the discretization of the multi-media mathematical model was completed: Oil phase: Water phase: Gas phase: In the formula, both the change of conductivity ΔT n+1 and the change of pressure Δp n+1 are the solution variables at n + 1 time. There are two terms of product at the left end of the equation, which is nonlinear. In order to realize the solvability of the equation, the equation needs to be linearized and expanded. When the fully implicit method is used to solve the equations, each time step needs to go through a series of iterations.

Quick Solution Method
At present, the calculation method embedded in threedimensional three-phase numerical simulation software is usually the ESPIDO algorithm. When dealing with the edge matrix, this method adopts the method of matrix conversion of the edge matrix and eliminates the edge elements to solve the problems of solving the edge matrix. The specific matrix elimination method is shown in Figure 5.
It can be seen from Figure 5 that this method of dealing with the edge matrix largely depends on whether the matrix in the lower right corner can be converted into a unit matrix. This transformation form will not bring much trouble when dealing with the traditional edging problems (such as implicit bottom hole pressure problem and double hole single permeability problem). However, it is prone to the problem of nonconvergence or small convergence radius when facing the problem of multiple medium coupling. By investigating the commonly used large sparse matrix solution methods (Jianhua, 2007;Jianhua and Hua, 2008), the advantages and disadvantages of different algorithms were determined, and the adaptability of different algorithms to the solution of multi-media coupling model of tight reservoirs was evaluated, as shown in Table 1.
Considering the calculation workload, convergence stability, residual norm, and convergence speed, the generalized Gl-CGS algorithm with strong applicability was the generalized global conjugate gradient square algorithm (Zhongxiao, 1998).

APPLICATION OF THE MATHEMATICAL MODEL
Herein, new numerical simulation software for tight reservoirs was applied to study the tight reservoirs in Wang 587 Block of NiuZhuang oilfield in Jiyang Depression of Bohai Bay Basin. Firstly, the calculation accuracies of the development index of the new software and the commercial Eclipse software were compared. Then, the reasonable development policy was optimized for the problems of low liquid production and poor reservoir production effect in the Wang 587 block, which provided a theoretical basis and scheme guidance for the development scheme design of tight reservoirs in turbidites in other blocks.
Based on the seismic data, logging data, and laying data of Wang 587 Block, fine geological modeling was carried out, including the structural model, sedimentary facies model, and property model. The top structure diagram of the No. 3 sand formation, member 4 of Shahejie Formation in Wang 587 Block, is shown in Figure 6.
The structural model: the maximum principal stress direction in this block is NE72°. In order to facilitate the design of compression crack and avoid grid effects, the grid was divided along the principal stress direction with a plane grid step of 20 × 20 m. Vertically, taking the single sand body as the unit, the internal subdivision grid step of the single sand body was 0.2 m, the longitudinal heterogeneity was simulated, and the structure model was established by well control thickness modeling, as shown in Figure 7. The sedimentary facies model: the sedimentary phase can truly reflect the river dynamics and the sedimentary rhythm. Its spreading mode has an important connection with the property distribution, which is an important means to accurately control the non-homogeneous characteristics of various sand bodies. The sedimentary facies model is the basis for establishing the facies-controlled property model. After importing the discrete sedimentary facies data, the sedimentary facies points were processed into sedimentary facies surfaces, and the sedimentary facies surfaces were integrated to establish the sedimentary facies model, as shown in Figure 8. Facies controlled property model: Wang 587 Block has dense well points and complete logging and layering data, suitable for the deterministic modeling method. The property model of Wang 587 Block was established by using the facies-controlled property modeling method, including the permeability model and oil saturation model, as shown in Figure 9. The property model was calculated based on interpolation of the same sedimentary facies type, considering the superior and inferior characteristics of the sedimentary facies and randomly interpolating within a reasonable range, which preserved the randomness and ensured the reasonableness of the properties.
Compared with the Eclipse model, the applicability of the nonlinear seepage coupling model to tight reservoirs was mainly reflected in the calculation accuracy of oil production in the initial stage of development. Herein, the Wang 587 Block numerical simulation model was established using different software. Based on the actual development data of the single well, the accuracy of the calculation results of the new model and Eclipse numerical simulation was verified. The reasons for the difference in calculation results were analyzed. The Petrel property model was coarsened and imported into the structural module of the numerical simulation. The numerical simulation model of Wang 587 Block was established using the new numerical simulation software and commercial software Eclipse. Two working systems were mainly used for numerical simulation: constant liquid production and constant pressure production. During constant liquid production, if the liquid production was lower than the set liquid production, the bottom hole flowing pressure in the software would automatically decrease until the standard atmospheric pressure was 0.1 MPa. Thus, there would be the phenomenon that the calculated liquid production was consistent with the actual flow pressure, but the bottom hole flowing pressure was lower than the actual flow pressure. Moreover, there would be the phenomenon that the calculated liquid production was lower than the actual liquid production, and the bottom hole flowing pressure was reduced to 0.1 MPa. The accuracy of bottom hole flowing pressure can be ensured during the production of constant bottom hole flowing pressure. On this basis, the calculated and actual liquid production can be compared to achieve the purpose   of controlling variables. Therefore, the working system of constant bottom hole flowing pressure in this article was adopted to predict the development indicators. Taking the actual initial development data of well W587-X3 as an example, the calculation effects of liquid production, water cut, and oil production were compared, as shown in Figures  10-12.
The oil production validation results of different models showed that ① the calculation accuracy of the nonlinear seepage coupling model was 95.6%, which was in high agreement with the reality. While the average fitting rate of the Eclipse model was only 82.3%, especially in the early elastic development stage, the calculation error was as high as 34.1%. ② The low initial fluid production results and high late fluid production results calculated in the Eclipse software were because oil production from fracturing development in tight reservoirs was mainly influenced by the matrix/microfracture transfer flow, as shown in Figure 13. The transfer flow quantity of the transfer flow model used in the nonlinear seepage coupling was large at the initial stage but decreased rapidly. After entering the pseudo-steady flowing period, the transfer flow was small, which was highly consistent with the actual transfer flow. The calculation results of the steady-state transfer flow model used in Eclipse were lower in the early stage and higher in the late stage because the influence of elasticity was not considered in the early stage and the liquid production was low. The influence of threshold pressure gradient was not considered in the late stage, which exaggerated the production pressure difference and high liquid production. Generally speaking, the Eclipse model cannot simulate the phenomenon of large production decline in the early stage of tight reservoir development, and the calculation error of liquid production was large after entering the quasi-stable period.
The validation results of different models for fluid production and water cut showed that the nonlinear seepage coupling was consistent with the actual tight reservoir development state due to the consideration of the threshold pressure gradient, pressure-sensitive effect, and matrix/microfracture unsteady transfer flow. Moreover, the application of nonlinear seepage coupling numerical simulation software to optimize the development scheme was highly credible.

CONCLUSION
(1) The multi-media coupling seepage mathematical model of matrix/microfracture/main fracture for fracturing development of tight reservoirs was established. The effects of nonlinear seepage, threshold pressure gradient, pressure-sensitive effect, and matrix/fracture nonlinear transfer flow were considered, which can fully reflect the fracturing development characteristics of tight reservoirs.   (2) Based on the fully implicit simultaneous solution, the multimedia coupling discrete mathematical model was constructed, and the fast solution algorithm of the generalized global conjugate gradient was optimized, which provided theoretical support for the development of tight reservoir numerical simulation software.
(3) The model application effect of tight reservoirs in Block Wang 587 showed that the fitting accuracy of the new tight reservoir numerical simulation software was as high as 95.6%, which was 13.3% higher than that of Eclipse software, and the reliability of using the new numerical simulation software to optimize the tight reservoir development scheme was higher.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding author.