A New Approach for Automatic Optimization of Complex Fracture Network in Shale Reservoirs

Shale gas reservoirs have gradually become the main source for oil and gas production. The automatic optimization technology of complex fracture network in fractured horizontal wells is the key technology to realize the efficient development of shale gas reservoirs. In this paper, based on the flow model of shale gas reservoirs, the porosity/permeability of the matrix system and natural fracture system is characterized. The fracture network morphology is finely characterized by the fracture network expansion calculation method, and the flow model was proposed and solved. On this basis, the influence of matrix permeability, matrix porosity, fracture permeability, fracture porosity, and fracture length on the production of shale gas reservoirs is studied. The optimal design of fracture length and fracture location was carried, and the automatic optimization method of complex fracture network parameters based on simultaneous perturbation stochastic approximation (SPSA) was proposed. The method was applied in a shale gas reservoir, and the results showed that the proposed automatic optimization method of the complex fracture network in shale gas reservoirs can automatically optimize the parameters such as fracture location and fracture length and obtain the optimal fracture network distribution matching with geological conditions.


Introduction
Conventional natural gas resources have gradually dried up after long-term exploitation, and unconventional oil and gas resources such as tight oil and gas, shale gas, coalbed methane, and gas hydrate have gradually become the main source for increasing oil and gas reserves and production [1]. Shale gas is an important part of unconventional oil and gas resources. The development practice of shale gas reservoirs in many areas has proved that fractured horizontal wells are the main technology for effective exploitation of shale gas reservoirs [2]. Accurate and efficient acquisition of optimal fracture parameters can provide help for oil and gas field decision makers.
The stimulated shale gas reservoirs mainly include matrix, natural fractures, and hydraulic fractures. Shale matrix is dense; pore size is very small, which distributes organic matter and inorganic matter; and fluid occurrence form is complex. Shale gas mainly exists in the surface of organic matter as an adsorption state and in the intergranular pores as an adsorption state and free state. The nanoscale pores in shale lead to a large internal surface area, and the content of adsorbed gas may exceed 50% [3]. There is bound water film on the internal surface of inorganic pores, and shale gas is mainly in a free state, followed by an adsorption state. Scholars have proposed many models to describe the apparent permeability of shale gas reservoirs [4]. Javadpour [5] divides the microscopic flow into two parts, namely, the flow caused by Knudsen diffusion and the traditional flow under pressure difference. It is considered that the total flow is equal to the flow caused by Knudsen diffusion plus the flow caused by pressure difference. Then, the flow formula is compared with the flow formula obtained by Darcy law, and a parameter like the permeability is deduced, that is, the apparent permeability. By introducing the apparent permeability, Darcy law is still applicable to the micropore flow. Sakhaee-Pour and Bryant [6] introduced the slip boundary to correct the apparent permeability, and the correction formula was given based on the Beskok-Karniadakis slip boundary. The simulation results are verified by the data of experiment, which proves that the new model can well match the experimental data.
Natural fractures and hydraulic fractures coexist in shale gas reservoirs, where natural fractures are randomly distributed and hydraulic fractures have complex distribution patterns, which are difficult to be accurately characterized. Conventional fracture characterization methods often assume that the fractures in the reservoirs are evenly distributed [7], which is very different from the fracture distribution in the actual shale gas reservoirs, resulting in this kind of model cannot accurately describe the fracture distribution in shale gas reservoirs. According to various exploration and production data, the distribution of fractures in shale gas reservoirs is very uneven, and the porosity and permeability in shale matrix are very low. Therefore, fractures in shale gas reservoirs are the main channel for gas flow, and the distribution and properties of fractures have a great impact on the productivity of shale gas reservoirs [8]. Consequently, the numerical simulation for shale gas reservoirs must overcome the problem that the conventional methods cannot accurately describe fractures. In 2020, Zhao et al. [9] used the similarity principle to compare the similarity between the lightning formation process and the fracture network expansion and firstly proposed a calculation method for the fracture network expansion of reservoir fracturing based on lightning simulation. This method can quickly conduct large-scale fracture propagation simulation, which provides the possibility for the fracture propagation simulation of the entire reservoirs.
Multiscale porous media in shale gas reservoirs have strong heterogeneity, and the fluid migration mechanism of different scale media is complex. It is of great significance to clarify the multimedia fluid cross-scale mass transfer characterization method for shale gas reservoir numerical simulation. At present, the research methods of the multiscale medium coupling flow model are mainly divided into continuous medium model and discrete fracture model [10][11][12][13][14][15]. Noorishad and Mehran [16] and Baca et al. [17] proposed the discrete fracture model (DFM model) for the first time, ignoring the change of the physical quantity of the fracture in the opening direction. The mesh in the fracture is combined by the geometric method and solved by the finite element method. For N-dimensional reservoirs, fractures can be treated as N − 1 dimensional. In a three-dimensional reservoir system, fractures are treated as a two-dimensional plane, while in a twodimensional system, fractures are considered as a onedimensional line segment. This process can significantly reduce the amount of computation and computation time [18][19][20][21]. Kim and Deo [22] established a discrete fracture multiphase flow model. The fractures in the model are randomly distributed and can intersect, which is applied to two-dimensional water flooding simulation. The matrix is divided by a triangular mesh, and the fracture is divided by a linear element segment mesh, which is solved by the Galerkin finite element method. Wu et al. [23] proposed a fractal discrete fracture model. This method is based on fractal theory, by changing the fractal parameters to control the probability distribution of fracture center, length, and density.
After the establishment of a numerical simulation method for shale gas reservoirs, it is necessary to optimize the perforation parameters of horizontal wells for specific reservoir geological conditions. The general optimization method mainly uses artificial scheme design combined with numerical simulation to optimize the fracturing design [24]. Bu [25] adopted the artificial optimization method and selected the highest production by numerical simulation method. In the Xinchang gas field, the nonuniform fracture length, nonuniform fracture spacing, fracturing scale, and fracture number are optimized. Cao et al. [26] established the CO 2 flooding geological model of the inverted sevenpoint horizontal well pattern and optimized the parameters of the inverted seven-point horizontal well pattern. Zhang and Jiang [27] analyzed the influence of injection production well spacing, horizontal section length, fracture half-length, and other parameters and obtained the optimal combination of fracture network parameters. Artificial scheme combined with the numerical simulation method cannot accurately obtain the global optimal parameters, well pattern and fracture pattern parameters cannot be optimized synchronously, and artificial scheme design is time-consuming and labor-consuming, which could not realize the realtime optimization for reservoir development. Huang et al. [28] introduced intelligent optimization technology into the application of horizontal well fracturing development and obtained more accurate and flexible results than orthogonal test optimization, which caused extensive research by scholars. Xu et al. [29] proposed a framework of the embedding discrete fracture model (EDFM) and intelligent algorithm to optimize the parameters of fractured horizontal wells and applied it in actual reservoirs to obtain a higher net present value than local grid refinement (LGR). However, previous studies mainly focused on the optimization of fracture parameters in the fixed fracture pattern [30][31][32][33][34][35][36], which could not automatically optimize the fracture pattern distribution of fractured horizontal wells according to geological conditions.
Many scholars have carried out integrated optimization research on well pattern and fracture network, so that the results are more in line with the understanding of geological characteristics. Feng et al. [37] took the maximization of economic net present value as the optimization objective; first optimized the three fracture network parameters of conductivity, spacing, and half fracture length; and then established three operators to optimize the well network parameters, which improved the optimization efficiency without reducing the accuracy of the optimization results. Ma et al. [38] proposed an integrated optimization method of well pattern and fracture parameters based on the SPSA algorithm and realized automatic continuous optimization for key parameters such as horizontal well position, horizontal well length, and fracture parameters. An example was applied to obtain the optimal parameters of well pattern and fracture parameters matching with geological conditions. However, the existing integrated optimization methods are not detailed enough in some specific problems, such as the inaccurate flow model, the optimization parameters not taking into account all the conditions, and the unclear characterization of the fracture network shape.
In view of the existing problems, this paper fully considers the characteristics of flow model for shale gas reservoirs. Based on the latest fracture propagation simulation method, a flow simulation method considering the morphological characteristics of complex fracture networks is established and solved by the finite volume method. On this basis, the SPSA algorithm is used to automatically optimize the hydraulic fracture length and fracture location.

Flow Model for Complex Fracture
Network in Shale Gas Reservoirs  [39], the microunits in the reconstruction area are intercepted, and the physical flow model of orthogonal discrete fracture coupled with dual-media matrix block is established [40]. As shown in Figure 1(b), the fracture network distribution in the actual reservoir reconstruction area is relatively complex, and it is typical as a vertical cross fracture network. The fracture network divides the reservoirs into cube matrix blocks. The physical properties of matrix blocks are the same as the initial properties of shale gas reservoirs; that is, organic matter and inorganic matter are developed in matrix blocks (Figure 1(c)). Inorganic matter mainly contains free phase gas molecules. Micro-and nanopore throats are developed in inorganic matter. Gas migration in inorganic matter mainly includes Knudsen diffusion and viscous flow. Nanoscale organic intragranular pores are developed in the organic matter, and shale gas is mainly stored in the organic matter in the form of adsorption phase and free phase (Figure 1(d)). The gas migration in the organic matter mainly includes the diffusion on the surface of the organic matter, the Knudsen diffusion, and viscous flow in the nanopores. Without considering the influence of in situ stress and other factors on fluid migration, the gas stored in organic matter will flow to the bottom of the well through the following paths: internal migration of organic matter, cross-scale mass transfer between organic matter and inorganic matter, internal migration of inorganic matter, cross-scale mass transfer between inorganic matter and fracture network, and internal migration of fracture network.

Characterization of Porosity/Permeability of Matrix
System. The matrix region assumes the discrete distribution of organic and inorganic matter, as shown in Figure 2.
The gas in organic matter occurs in the form of free gas and adsorbed gas. The gas migration in organic 3 Lithosphere includes the following forms: viscous flow considering slip effect, Knudsen diffusion, surface diffusion, and Langmuir adsorption/desorption. According to the previous research results, the apparent porosity/permeability can be expressed as [41] where k app k is the apparent permeability of organic matter coupling multiple migration mechanisms, mD; ϕ app k is the apparent porosity of organic matter coupling free gas and adsorbed gas; ϕ k is the inherent porosity of organic matter; p k is the pore pressure of organic matter, MPa; Z is the gas compression factor; R is a general gas constant of 8:314 × 10 −6 MJK −1 mol −1 ; T is the formation temperature, K; ε ks is the kerogen solid volume in unit core volume; M is shale gas molecular weight, kg/mol; μ g is gas viscosity, MPa s; r k is nanoscale pore diameter in kerogen, m; τ k is the pore tortuosity correction coefficient in kerogen; C g is the gas compression coefficient, MPa −1 ; D s is the diffusion coefficient of adsorbed gas on kerogen surface, m 2 /s; c μs is the maximum adsorption concentration of Langmuir on the surface of kerogen, mol/m 3 ; and p L is the Langmuir pressure, MPa. f ks is the kerogen angular momentum adjustment coefficient, which is related to the diffusion and reflection of shale gas on the nanopore wall in organic matter, and the variation range is 0-1.
The shale gas in inorganic matter mainly occurs in the form of free gas. Gas migration in inorganic matter micropore media includes viscous flow and Knudsen diffusion considering slip effect. Similarly, the apparent permeability    Lithosphere of inorganic matter coupled with gas viscous slip flow and Knudsen diffusion can be obtained [41]: k app m is the apparent permeability of inorganic matter coupled with multiple migration mechanisms, mD; φ m is the porosity of inorganic matter; p m is the pore pressure of inorganic matter, MPa; r m is the pore diameter of micron scale in inorganic matter, m; τ m is the correction factor of pore tortuosity in inorganic matter. f ms is the adjustment coefficient of inorganic matter angular momentum, which is related to the diffusion and reflection of shale gas on the micropore wall of inorganic matter, and the variation range is 0-1.
In the dual medium model, the discontinuous "matrix" and "fracture" in the discrete medium model need to be equivalent to the dual continuous medium. According to the basic assumption of the model, different medium size, and quantity distribution characteristics, the discrete medium model in Figure 2 is treated as continuity equivalent. The apparent porosity and permeability of the matrix are as follows [41]: k equ k is the equivalent permeability of organic matter in dual media, mD; f k is the characteristic size of organic matter in the discrete media model, m; n k is the number of organic matter clumps in the discrete media model; L KM is the side length of inorganic matter characterization unit model, m; φ equ k is the equivalent porosity of organic matter in the dual media; k equ m is the equivalent permeability of inorganic matter in dual media, mD; f m is the characteristic size of inorganic matter mass in the discrete media model, m; n m is the number of inorganic matter masses in the dis-crete media model; φ equ m is the equivalent porosity of inorganic matter in the dual media model.

Porosity/Permeability Characterization of Natural
Fractures. According to the characteristics of natural fractures in shale gas reservoirs, using tensor theory and equivalent flow resistance principle, the reservoir is composed of many fractures developed regions and matrix regions without fractures (Figure 3).
Assuming that the average width (opening) of natural fracture is b f , the linear density of fracture is D f , the total width of unit is b t , and the total width of rock fracture is b c = D f b t b f , the equivalent permeability coefficient of natural fracture system can be defined as 2.1.4. Characterization of Secondary Fractures. In 2020, Zhao et al. [43] proposed a calculation method of fracture network expansion based on lightning simulation. This method can describe the fracture network morphology more precisely and has the characteristics of high computational efficiency and large-scale fracture network inversion. The inverted fracture network morphology is shown in Figure 4. Based on the above fracture network morphology, the distribution morphology of fracture network can be accurately described. However, in the process of flow simulation, it is necessary to further determine the distribution of permeability and porosity of fractures to conduct flow simulation. The fractal diffusion equation proposed by Sheng et al. [44] is used to characterize the distribution of fracture porosity/permeability, namely, where k f is the equivalent permeability of fracture system, mD; k w is the fracture permeability at the reference point, mD; x w is the position of reference point in cm; x is the position of current point in cm; d fa is the fractal dimension of 5 Lithosphere induced-fracture aperture; ϕ f is the fractal porosity; and ϕ w is the porosity of fracture system at the reference point.

Model
Solving. Based on the specified fracture morphology and matrix porosity/permeability equation, the flow simulation is realized by a numerical method. Three kinds of media are developed in the reservoirs, matrix-natural fractures are characterized by dual media, and large-scale fractures are simulated based on an embedded discrete fracture model (EDFM). The equations of gas-water two-phase flow are discretized by a finite volume method, and the discrete equations are calculated efficiently and accurately by Newton iteration and automatic differential technology.

Conceptual Example.
In order to verify the correctness of the above mathematical model, based on the fracture propagation simulation method, this paper generates a complex hydraulic fracture network according to the fracture length and uses the fractal diffusion equation to characterize the distribution of porosity/permeability of hydraulic fractures. The flow simulation and solution are carried out for the matrix-natural fracture-hydraulic fracture system of shale gas reservoirs, respectively. The related reservoir parameters, fluid parameters, and fracture network parameters of the reservoirs are shown in Table 1. The specific   Lithosphere fracture morphology, pressure distribution, and production curve are calculated, as shown in Figure 5. WGPT in the figure is cumulative gas production.

Sensitivity Analysis.
Based on the above model, this paper explores the effects of different matrix permeability, matrix porosity, fracture permeability, fracture porosity, and fracture length on pressure distribution and production curve.
(1) Effect of Matrix Permeability. The matrix permeability will significantly affect the fluid seepage capacity. This paper selects three different matrix permeabilities, which are 0.01 md, 0.02 md, and 0.03 md, respectively, to explore the influence of matrix permeability on gas production, as shown in Figures 6 and 7. Compared with the calculation results, it is found that under the same conditions, with the increase of matrix permeability, the fluid flow ability from matrix to fracture network increases, and the shale gas production increases gradually. However, when the permeability is high, the formation fluid is easier to flow directly into the bottom hole, the control area of single fracture is closer to the round, and the fracturing effect is less obvious. This shows that in shale gas reservoirs with low matrix permeability, the effect of fracture network is better.
(2) Effect of Matrix Porosity. The size of matrix porosity represents the gas reservoir pore space, which directly affects the gas reserves. This paper selects three different matrix porosities, namely, 5%, 8%, and 12%, to explore the influence of matrix porosity on gas production, as shown in Figures 8 and 9. Compared with the calculation results, it is found that if the porosity of the matrix is higher, the geological reserves will be more abundant, and the shale gas production will gradually increase.
(3) Effect of Fracture Permeability. Fracture permeability refers to the fracture flow conductivity. If the fracture permeability is higher, the fluid in the fracturing area will flow into the wellbore more easily, but higher conductivity means more fracturing cost. This paper selects three different fracture permeabilities, namely, 10 mD, 20 mD, and 40 mD, to explore the influence of fracture permeability on gas production, as shown in Figures 10 and 11. It can be seen that with the increase of fracture permeability, the gas production in the fracture area increases gradually, while the gas inflow into the fracture in the formation basically remains unchanged. The comparison calculation results show that the fracture permeability has little effect on the gas production. When the fracture permeability is low, it is obvious that the fractal diffusion equation is more suitable for the actual situation. That is to say, the closer the wellbore is, the better the fracturing effect is, and the farther away the wellbore is, the worse the fracturing effect is, which realizes the effect of fine characterization of fracture permeability and porosity distribution.
(4) Effect of Fracture Porosity. Fracture porosity refers to the size of diversion space of fracture. If the fracture porosity is larger, the fluid in the formation is easier to flow into the fracture stimulated area, but higher porosity also means more fracturing cost. In this paper, three different fracture porosities, 1%, 8%, and 20%, are selected to explore the influence of fracture porosity on gas production, as shown in Figures 12 and 13. Compared with the 7 Lithosphere calculation results, it is found that the larger the porosity of fracture reconstruction area is, the greater the gas production is. The change of fracture porosity has a certain influence on the gas production in the fracture reconstruction area, but the effect is not obvious in the lowpermeability matrix area.
(5) Effect of Fracture Length. Fracture length is closely related to production. If the fracture length is longer, the effective fracturing volume of the well will be larger, but the long fracture length means high fracturing cost. In this paper, three different fracture lengths, 100 m, 140 m, and 180 m, are selected to explore the impact of fracture length on gas production, as shown in Figures 14 and 15. The results show that the higher the fracture length is, the larger the reconstruction area is, and the shale gas reservoirs with longterm and large-scale development will have higher production.

Optimization of Fracture Parameters Based on SPSA
3.1. Objective Function and Optimization Variables. In order to reduce input and increase output, this paper selects net present value (NPV) as the optimization objective function in shale gas reservoirs. Under different assumptions, NPV has different mathematical expressions, but it is usually composed of two parts [45], which are discounted cash flow and capital expenditure. The calculation formula of the objective function is as follows: In the formula, NPV is the economic net present value, yuan. N t is the total simulation time, day. r gas is the gas price, yuan/m 3 . Q gas ðtnÞ is the cumulative gas production in t n period, m 3 . b is the annual interest rate of banks, decimal. FC is the fixed cost, yuan/mouth. C well is the cost of horizontal well drilling, yuan/well. N f is the total number of fractures. C f i is the fracture cost of fracture i, yuan. According to the previous sensitivity analysis, when the reservoir geological characteristics are known, the number and location of fractures have the greatest impact on the gas production. Therefore, the length and location of fractures are used as the optimization variables for this optimization.

Principle of Simultaneous Perturbation Stochastic
Approximation. Simultaneous perturbation stochastic approximation (SPSA) obtains the search direction through synchronous disturbance of all control variables, and the direction is always the upward direction [46]. The calculation process only needs to solve the objective function value twice. Compared with the gradient method, the calculation amount is greatly reduced, which is very suitable for the optimization of high-dimensional problems. The basic principle of the SPSA algorithm is as follows.
A symmetric Bernoulli distribution vector δ k with parameter ±1 is generated. The probability of +1 or -1 is 50%, and the expected value is 0, so δ k −1 = δ k . In the k 8 Lithosphere In the formula, u k is the vector of the optimal control variable corresponding to step k. ε k is the disturbance step length. Jðu k Þ is the objective function value when the control variable is u k .
After obtaining the disturbance gradient, new control variables can be generated: In the formula, λ k is the search step of step k, and k·k ∞ is the infinite norm.

Optimization Process.
On the basis of given reservoir porosity and permeability, the fracture network expansion method is used to describe the fracture network morphology of different lengths, and the fractal diffusion equation is used to characterize the distribution of fracture porosity/permeability, and the numerical calculation method is used to realize the flow simulation solution. Then, the number of hydraulic fractures is set, the length and position of each hydraulic fracture are initialized, and the model output under the optimized variable is obtained by solving the model, and the corresponding objective function value is calculated. Further, the SPSA optimization algorithm is used to automatically optimize and adjust the optimization variables until the objective function converges or reaches the maximum number of iterations. The specific optimization process includes the following four steps: (1) Setting the initial fracture length and fracture position of each fracture, the gas reservoir numerical simulation model is established, and the initial gas production and initial objective function value NPV0 are obtained (2) According to the initial variable value, the SPSA algorithm is used to calculate the disturbance gradient from Equation (10), and the new control variable is obtained by updating the variable through Equation (11) (3) According to the updated variable values, the model file is automatically modified by the program, and the objective function value NPV is obtained by establishing and solving the new model 9 Lithosphere (4) If the NPV obtained in step 3 is always equal to the NPV0 of the previous iteration step, or the number of iterations exceeds the maximum number of iterations, then proceed to the next step. Otherwise, replace the initial variable and initial NPV0 with the updated variable and objective function NPV, and repeat steps 2-4 (5) Output optimal variables and objective function values

Application Example
In order to verify the correctness of the optimization method, a two-dimensional heterogeneous gas reservoir model is selected for demonstration ( Figure 16). The geological model parameters are as follows: the model size is 500 m × 400 m × 5 m, grid step is 1 m × 1 m × 1 m, the reservoir depth is 2000 m, the effective thickness is 5 m, and the average porosity is 0.08. The average original gas saturation is 20%, the original formation pressure is 25 MPa, and the surface gas density is 0.79 kg/m 3 . Horizontal wells are placed horizontally in the middle of the gas reservoir, the halflength of fractures is 100-180 m, and the number of hydraulic fractures is fixed to 6.
The relevant parameters of production cost are set as follows: the total production time is 10 years, the annual interest rate is 10%, the cost of hydraulic fracturing is 8 million yuan, the fixed cost is 1 million yuan, the perforation cost is 100000 yuan/piece, and the price of natural gas is 2 yuan/m 3 , regardless of the cost of water production.
The location and length of hydraulic fractures are optimized. When the optimization reaches the 36th step, the convergence is achieved. The changes of optimization variables and economic net present value of some iteration steps are shown in Table 2, where Loc 1~L oc 6 are the horizontal position of pressure fractures. The economic net present value curve shown in Figure 17 is made according to the economic net present value of each iteration step.
It can be seen that compared with the parameters of the initial fractured horizontal well, the final fractured horizontal well has changed greatly, and the economic net present value has been optimized from 7.0655 million yuan to 16.6410 million yuan, an increase of 9.5755 million yuan. The overall porosity of the reservoir is the same, the permeability is relatively high in the middle and upper parts, and the physical properties of the formation are better. Therefore, the final optimization result is that a shorter fracture is deployed in the middle part, which can obtain better production at a lower cost. However, the permeability of the lower left part is relatively low, and the gas is not easy to flow  10 Lithosphere into the wellbore, so the longer fractures are used for exploitation, so as to improve the natural gas production. The optimization results accord with the actual understanding of mineral resources, which proves that the method can be applied to the actual reservoirs.

Conclusions
(1) For the physical model of fractured horizontal wells, considering the permeability/porosity characterization of matrix and natural fractures, the fracture network morphology is finely described, and the Newton iteration and automatic differential technology are used to efficiently and accurately calculate the discrete equations (2) The sensitivity analysis shows that if the matrix permeability, matrix porosity, fracture permeability, fracture porosity, and fracture length are larger, the cumulative gas production will be more. In shale gas reservoirs with lower matrix permeability and higher porosity, the improvement effect of fracture network is more obvious. From the point of view of economic benefits, there is an optimal value between fracture conductivity and fracture length (3) It is proved that the automatic optimization method of complex fracture network in shale gas reservoirs can automatically optimize the length and location of main fractures according to the different geological characteristics of gas reservoirs and obtain the fracture network distribution with the highest economic benefit

Data Availability
The data would be supplied on request.