Ensemble-based optimization of hydraulically fractured horizontal well placement in shale gas reservoir through Hough transform parameterization

Shale gas reservoirs have been successfully developed due to the advancement of the horizontal well drilling and multistage hydraulic fracturing techniques. However, the optimization design of the horizontal well drilling, hydraulic fracturing, and operational schedule is a challenging problem. An ensemble-based optimization method (EnOpt) is proposed here to optimize the design of the hydraulically fractured horizontal well in the shale gas reservoir. The objective is to maximize the net present value (NPV) which requires a simulation model to predict the cumulative shale gas production. To accurately describe the geometry of the hydraulic fractures, the embedded discrete fracture modeling method (EDFM) is used to construct the shale gas simulation model. The effects of gas absorption, Knudsen diffusion, natural and hydraulic fractures, and gas–water two phase flow are considered in the shale gas production system. To improve the parameter continuity and Gaussianity required by the EnOpt method, the Hough transformation parameterization is used to characterize the horizontal well. The results show that the proposed method can effectively optimize the design parameters of the hydraulically fractured horizontal well, and the NPV can be improved greatly after optimization so that the design parameters can approach to their optimal values.


Introduction
Shale gas development has received great attention in recent years because of the increasing demand for natural gas resources and its huge reserves worldwide (Dong et al. 2012(Dong et al. , 2014. The matrix permeability of the shale gas reservoir is ultra-low due to the nanoscale pore size (Javadpour et al., 2007;Liu et al. 2019a). With the advancement of horizontal well drilling and hydraulic fracturing techniques, the natural gas resources in shale reservoir can be successfully unlocked in a cost-effective manner.
Shale gas is considered as one of the unconventional energy resources, and many special physical mechanisms are associated with the gas flow in shale reservoirs. The natural gas in the shale gas reservoir can mainly exist in two forms, i.e., free gas and adsorbed gas (Curtis 2002). The adsorbed gas will be desorbed into the matrix pores from the organic matters when the reservoir pressure decreases during the production process. In the shale matrix, the non-Darcy flow can occur due to the nanoscale pore size, and the Knudsen diffusion has been used to account for the effect of the molecular scale gas flow (Javadpour 2009;Xue Edited by Yan-Hua Sun et al. 2019;). In the fracture system of the shale gas reservoir, the natural fractures and the hydraulic fractures together form a complex fracture network, and the shale gas will flow through these fractures with various sizes into the horizontal well (Cipolla et al. 2010;Clarkson et al. 2011). The multistage hydraulic fracturing will create multiple transverse fractures along the lateral direction of the horizontal well in relatively intact reservoirs, or generate complex fracture networks connected to the horizontal well in fractured reservoirs (Li et al. 2016Li and Zhang 2018;Liu et al. 2019b). To improve the ultimate recovery of the shale gas in a cost-effective manner, the design of the hydraulically fractured horizontal well needs to be optimized.
Fractures in the hydraulically fractured shale gas reservoir increase the contact surface between shale matrix and fractured horizontal well and form the primary flow path for shale gas flow. Reservoir numerical simulation and datadriven methods can forecast the shale gas production and aid the optimization of the fractured horizontal well Xue et al. 2020Xue et al. , 2021. To take the fractured flow into consideration, multi-continuum models, such as the dual porosity model (Warren and Root 1963;Kazemi 1969;Rostami et al. 2020) and multi-porosity model (Wu and Pruess 1988), have been used for shale gas reservoir simulation. Zhang et al. (2009) conducted the sensitivity analysis of the hydraulic fracture parameters in the shale gas reservoir by upscaling the discrete natural fractures to dual porosity system. Cipolla et al. (2010) simulated the shale gas production by using the dual porosity model with the consideration of gas adsorption effect. Rubin (2010) used local grid refinement (LGR) method with dual porosity model to improve the characterization of the hydraulic fractures. Essentially, the multi-continuum model represents the fractures by using an averaged property in the grid block. The discrete fracture modeling method can characterize the fracture geometry explicitly (Karimi-Fard et al. 2004;Gong et al. 2008), but the unstructured grid has to be used and the computational burden can be too heavy to be implemented in the fieldscale application. The embedded discrete fracture modeling method (EDFM) was proposed by Li and Lee (2008) and improved by Moinfar et al. (2014), which can represent the geometry of each fracture explicitly but with structure grid system, and this method has been applied to simulate the shale gas production (Dai et al. 2017;Yu et al. 2018). Due to the uncertainty associated with the geological condition including the conceptual model uncertainty (Xue and Zhang 2014;Xue et al. 2015), the production data can be used to reduce the predictive uncertainty and obtain more accurate prediction results (Dachanuwattana et al. 2018;Kang et al. 2020;Xue et al. 2020).
Many optimization methods have been investigated to deal with the well placement problem in reservoir engineering. Brouwer and Jansen (2004) used a systematic dynamic optimization approach to control the valve setting in the injector and producers for a water flooding project in the heterogeneous reservoir. Sarma and Chen (2008) and Zhang et al. (2010) converted the discrete control variables to continuous ones and used the adjoints and gradientbased optimization method to determine the optimal the well locations. Volkov and Bellout (2018) approximated the well placement optimization gradient by using finite difference approximations of augmented Lagrangian derivatives with adjoint formulation to improve efficiency. Al Dossary and Nasrabadi (2016) proposed an imperialist competitive algorithm (ICA) to optimize the well location. Hamida et al. (2017) modified the traditional genetic algorithm with a similarity operator to optimize the well placement. Jesmani et al. (2015) used the derivative-free particle swarm optimization algorithm to optimize the well location with well design constrains. Tukur et al. (2019) implemented the genetic algorithm and simulated annealing methods to optimize the well placement. In the optimization of a hydraulically fractured horizontal well in a shale reservoir, Ma et al. (2013) optimized the hydraulic fracture placement with the gradient-based finite difference method (FD), discrete simultaneous perturbation stochastic approximation (DSPSA), and genetic algorithm (GA). Yu and Sepehrnoori (2013) used the response surface method to conduct the optimization of the multiple hydraulically fractured horizontal wells in shale gas reservoirs. Wilson and Durlofsky (2013) introduced a direct search algorithm to optimize the shale gas field development by combining with reduced physics model to improve computational efficiency. Rammay and Awotunde (2016) used differential evolution algorithm to optimize the hydraulic fracturing parameters and horizontal well length. Zhang and Sheng (2020) optimized the fractured horizontal well in the shale gas reservoir by considering the effect the stimulated reservoir volume. The ensemble optimization method (EnOpt) is powerful optimization method proposed by Chen et al. (2009). It is a stochastic gradient optimization method, which approximates the gradient through the ensemble computation, and thus it is capable of integrating with any simulator. In addition, compared to other gradientbased optimization method (e.g., finite difference method and adjoint method) and gradient-free optimization method (e.g., genetic algorithm and particle swarm method), it is more efficient when the dimension of the control variables is high. Leeuwenburgh et al. (2010) applied this method to optimize the settings of inflow control valve in a water flooding project and compared its performance with the adjoint method. However, this method has not been fully investigated in the optimization problem of hydraulically fractured horizontal well in shale gas reservoir.
In this research, we propose to use ensemble optimization method to maximize the net present value (NPV) in the shale gas development project. The EDFM is used to simulate the shale gas production. The optimization parameters are the location of horizontal well drilled in a shale reservoir with considering the effect of the natural fractures, the optimal hydraulic fracturing parameters (i.e., number of stages, length of the hydraulic fractures, and conductivity of hydraulic fractures) and the operational parameter (i.e., bottom hole pressure). Usually, these parameters are optimized independently step by step. However, in our research, all these parameters are optimized simultaneously so that the most optimal shale gas development plan can be designed. When optimizing the horizontal well location, the coordinates of the well are discrete variables. Here, we further propose to use Hough transformation-based parameterization to transform the Cartesian space into Hough space. This transformation can improve the Gaussianity of the parameter so that it can be better integrated with EnOpt method.

Shale gas flow simulation
The EDFM is used to solve the simulation of shale gas production. The detailed equations to characterize the multiscale shale gas flow have been presented in Xue et al. (2020), and the design of the simulation software was introduced in Li et al. (2015). The brief governing equations of the gas-water two phase flow model of shale gas production through EDFM method are introduced here.
The gas flow equation in shale matrix is: where g is the density of the shale gas, k m a is the apparent permeability of shale matrix, k m rg is the relative permeability of gas in the shale matrix, g is the gas viscosity, m is the matrix porosity, S m g is the gas saturation in the shale matrix, q s c is the desorption/adsorption mass flow rate, which can be characterized by Langmuir isothermal adsorption model, and q mf g is the mass flow rate from matrix to fractures. To account for the Knudsen diffusion, the apparent permeability k m a can be evaluated by (Tang et al. 2005): where k m ∞ is the absolute permeability, C 1 and C 2 are constant coefficients, and Kn is the Knudsen number.
The shale gas flow equation in fractures is: where k f is the absolute permeability of fractures, k f rg is the relative permeability of shale gas in fractures, p f g is the gas pressure in fractures, q well g is the mass flow rate of shale gas to well, q mf g is the mass flow rate of shale gas from matrix to fractures, f is the fracture porosity, S f g is the gas saturation in fractures, and is a coefficient to account for non-Darcy flow.
The water flow equation in the shale matrix is: where w is the water density, k m is the apparent permeability of water in the shale matrix, k m rw is the relative permeability of water in the shale matrix, w is the water viscosity, p m w is the water pressure in the shale matrix, q mf w is the mass flow rate from matrix to fracture, and S m w is the water saturation. The water flow equation in the fracture is: where p f w is the pressure of water phase in fractures and q well w is the mass flow rate of water to well. In the EDFM, the mass flow rate from matrix to fractures, q mf g , and two fracture segments connections in the same matrix block or in two adjacent matrix blocks can be computed through the non-neighboring connection (NNC) (Moinfar et al. 2014): where p f g is the gas pressure in fractures, N nnc is the number of the non-neighboring connection grids, T nnc is the NNC transmissibility factor, A nnc is the contact area of fracture and matrix, k nnc is the harmonic mean of matrix and fracture permeability, and d nnc is the characteristics distance between matrix block and fracture plane.
Together with the NNC between matrix and facture in EDFM mentioned above, two more NNCs need to be considered, i.e., the two fracture segments connections in the same matrix block and two adjacent matrix blocks. The NNC transmissibility, T nnc , can be written as: where L int is the length of the fracture intersection, k f is the fracture permeability, w f is the fracture aperture, and d f is the normal distance between the center of the fracture and the fractures intersection. The developed EDFM model has been validated against the commercial software CMG to show its accuracy (Dai et al. 2017).

Ensemble optimization method
The ensemble optimization method was proposed by Chen et al. (2009), and it has been applied in several production optimization problems (Chen and Oliver 2010;Fonseca et al. 2014;Tueros et al. 2018). A vector of the control variable is defined in the EnOpt method, which contains all the variables that need to be optimized. The control variable vector can be defined as: where N x is the number of control variables.
The NPV is used as the objective function during the optimization process here, and the formula is written as: where i is the time step index, N t is the total production time, Q(i) is the cumulative shale gas production within the given time step which can be computed from the above EDFMbased shale gas production model, P g is the price of shale gas which is set as 3 CNY/m 3 , r is the discount rate which is set as 6%, w L is the length of horizontal well, P hw is the drilling cost per unit well length which is set as 30,000 CNY/m, hf xf is the half-length of the hydraulic fractures, hf stage is the stage number of the hydraulic fractures, and P hf is the well complete cost per unit fracture length, which is set as 20,000 CNY/m. All these parameter setting values are inferred from a laboratory report to describe the shale gas development in Zhaotong shale gas field in China. During the optimization process, the NPV is maximized by optimizing the control variables.
In each optimization step, the control variable is updated through where k+1 is the updated control variable, k is the control variable before updating, k is the coefficient determining the updating step size, xx is the prior covariance matrix of the control variables, and k is the sensitivity of NPV function g(x) to the control variables.
The prior covariance matrix of the control variables xx is: where N e is the ensemble size.
The control variable matrix can be expressed as: x ij is the ensemble mean of the control variables.
The product term xx T k can be approximated by where xg is the covariance between X and G. It can be expressed as: And the NPV vector G can be expressed as: In the EnOpt method, we can substitute Eq. (13) in Eq. (10) and use xx as the filtering matrix. Then the optimization parameters are updated by

Hough transform parameterization for horizontal well
Hough transformation was proposed by Hough (1962), which was used to detect lines. Duda and Hart (1972) extended Hough transformation to detect any arbitrary objects in the image analysis. Hough transform converts the parameters from the Cartesian coordinate space to the Hough space. Hough space is basically an accumulator space, and it uses the voting method in the accumulation space to find the local maximum value to conduct feature detection. The accumulator to transform x-y Cartesian coordinate system space into ρ-θ polar coordinate system of Hough space can be expressed as: where A( , ) is recording the how many sinusoidal curves actually pass through the ( , ) point in the Hough space, I(x, y) is the point in the Cartesian coordinate space. (⋅) is the Dirac delta function, and it can be defined as: As shown in Fig. 1, the straight line y = −x + 5 in the Cartesian coordinate system can be transformed into a point in the Hough space The accumulator in the Hough space provides the information on the total number of sinusoidal curves that pass through the point ( , ) , and all the grid blocks along the line segment, characterized by ( , ) in Hough space, will contribute to the accumulator. The Hough transform has the capability to transform any line to the parameter set that can take a better continuity than the discrete point in the Cartesian coordinate (Lu and Zhang 2015;Yao et al. 2018); therefore, it can be used to represent a line segment in the Cartesian coordinate space. By using this method, the horizontal well is not characterized by its endpoint coordinates x 1 , y 1 , z 1 and x 2 , y 2 , z 2 , but it is represented by a parameter set (as shown in Fig. 2): where is the vertical line distance between the origin to line segment in the x-y projection plane; is the angle between x-axis and the vertical line in the x-y projection plane; D is the distance between the vertical line and the center point of the line segment in the x-y projection plane; L is the length of the line segment in the x-y projection plane; is the angle between x-axis and the line segment in the x-z projection plane; and is the angle between z-axis and the line segment in the y-z projection plane; is the angle between y-axis and the line segment in the x-y projection plane.
The horizontal well drilled in the shale gas reservoir, which can be regarded as the line segment. When defining the horizontal well with the Hough transformation-based parameterization method, the horizontal well is within the x-y plane in the coordinate system, and = 0 and = 0 in this case. In addition, it can be seen from the triangular relationship that = . Therefore, the horizontal well can be represented by the parameter set in the optimization problem:

Equivalent permeability conversion of natural fractures
In the shale gas reservoir, there may exist a large set of natural fractures, which results in the shale gas reservoir with severe heterogeneous permeability distribution (Khanal and Weijermars 2019). The hydraulic fractures can be represented explicitly with its full geometrical properties in the EDFM method. However, the number of the natural fractures is very large, which is infeasible to obtain the information on each natural fracture. Even we can know the exact distribution of the natural fractures, the simulation model can be too time-consuming to be used in practice. Here, we establish a method to convert the discrete natural fractures to their equivalent permeability. Let us define the azimuth angle of the fracture is , the dip angle of the fracture is , and permeability parallel to the fracture is K (as shown in Fig. 3). The permeability tensor of each fractures can be computed by (Song et al. 2019): When a number of N fractures exist in a single grid block, the permeability tensor can be expressed by: With consideration of the shale matrix permeability, the final permeability tensor can be computed by:

Results and discussion
In this research, the EnOpt method based on Hough transform is used to optimize the economic benefits of shale gas reservoir produced by fractured horizontal wells. The location parameters of fractured horizontal wells are transformed into Hough space, and the EnOpt algorithm is used for the integrated optimization of the design parameters. The Hough transformation-based parameterization can be suitable to improve the continuity of the discrete design parameters, such as the location or central point coordinate of the horizontal well, and can be helpful to the linearization of the nonlinear parameter in the set optimization algorithm. The transformed parameters can improve the Gaussian assumption of the EnOpt process and thus provide a more reasonable optimization performance.

Synthetic model construction of the shale gas reservoir
To demonstrate the workflow and performance of the proposed EDFM-EnOpt method, a synthetic model of shale gas reservoir simulation is constructed here. In the synthetic model, a hydraulically fractured horizontal well is located in a 3D shale gas reservoir. The shale gas simulation model is constructed using EDFM method. During the simulation process, it considers the special gas flow mechanisms in the shale gas reservoir. The simulation model is constructed based on gas-water two-phase flow model. In the shale matrix flow, it considers the adsorption/desorption of the shale gas on the organic content in the shale reservoir and the Knudsen diffusion effect of the shale gas flow caused by the nanoscale pore structure. In the fracture flow, it considers the gas flow in natural fractures by upscaling the natural fractures using equivalent permeability method and the gas flow in hydraulic fractures by explicitly taking the fracture geometry into account.
The parameter values used in the synthetic model to characterize the shale gas reservoir properties are summarized in Table 1. These data are collected from a laboratory report to study the Zhaotong shale gas field in China and represent the state of the shale gas reservoir before the well drilling and gas production.
A large number of natural fractures can be developed in the shale gas reservoir. It is hard to obtain the fracture geometry for each natural fracture, but the statistical distribution of the natural fractures can be obtained by geophysical method, such as imaging well logging. Usually, the properties of the natural fractures follow a certain random distribution. For example, the natural fracture apertures follow a Gaussian distribution (characterized by the mean and variance in the generation function), the lengths of the natural fractures follow a uniform distribution (characterized by the minimal and maximal values in the generation function), and the azimuth follows a Gaussian distribution. The random discrete natural fracture model is generated by using the parameter setting listed in Table 2. The central point density of the natural fractures are set as 4.0 × 10 −8 / m 2 . According to the relationship of the generated fracture endpoints coordinate data and the computational grid, the fracture endpoint coordinate information is converted into computational grid property, i.e., the equivalent grid block permeability. Then, the shale gas reservoir geological model with natural fractures can be constructed as shown in Fig. 4.
For the hydraulic fractures, the number of hydraulic fractures is small and the fracture conductivity is large. They directly connect with the horizontal wellbore and the shale matrix, which form the primary flow path for shale gas and has a great impact on the flow field. Therefore, the hydraulic fractures should be accurately described by considering their explicit fracture geometry. In this paper, the hydraulic fractures are established by EDFM, and the hydraulic fractures are transverse fractures under the reservoir condition, that is, the hydraulic fractures are rectangular plates perpendicular to the horizontal stratum. The parameter values are listed in Table 3, which are used as the initial design of the hydraulic fracturing project before optimization.
Once combining the properties of shale matrix, natural fractures and hydraulic fractures, the reservoir model for the shale gas production simulation can be generated, as shown in Fig. 5.

Optimization results of the hydraulically fractured horizontal well
When the shale gas production simulation model is established, the EnOpt proposed in this research can be used to optimize the hydraulic fracturing parameters. The objective is to maximize the NPV shown in Eq. (9) by optimizing the parameters of well drilling, well completion and reservoir operator simultaneously. Traditionally, the optimization process is conducted in a sequential way where these parameters  are optimized one by one independently. This method cannot guarantee to find the optimal design scheme. Due to the ensemble feature of our proposed method and the joint optimization paradigm, the obtained solution is expected to be optimal globally. Before starting the optimization process, we need to prepare the parameters that will be optimized. The parameters to characterize the horizontal well placement are the parameters ( , , D, L) shown in Eq. (20). The well location is a crucial parameter that controls the well performance of the shale gas development. The EnOpt method requires a Gaussian updating during the optimization process. The initial realizations can be sampled from any prescribed distribution to avoid Gaussian assumption, but the updating process only uses the first moment (i.e., mean) and second moment (i.e., covariance) in the EnOpt method. If the distribution of the horizontal well parameter is far from the Gaussian distribution, the updating or the optimizing performance will be deteriorated step by step. In the EDFM simulation method, the horizontal well is characterized by each discrete grid block in the background matrix grid system. The Hough transformation-based parameterization method here is used to transform the discrete horizontal well grid block to the continuous Hough parameter space. The location of the horizontal well can be defined by the angle and radius (i.e., and ) in the Hough space. The parameters after Hough transformation distribution are more continuous and smoother, and thus the Hough transformation-based parameterization is beneficial to improve the performance of the EnOpt when they are combined together.
During the optimization process, the parameters to characterize the horizontal well with Hough transformationbased parameterization are updated through the optimization method. Once the horizontal well parameter updating is done and a shale gas simulation process is required to evaluate the NPV, the horizontal well parameters in the Hough space are transformed back to its original space by:  With these two endpoints, the horizontal well can be constructed into the matrix grid system by setting the grid blocks that intersect with the straight line between these two endpoints as horizontal well. The parameters to characterize the hydraulic fracturing design are the number of fracturing stage hf stage , the halflength of the hydraulic fractures hf xf , and the conductivity of hydraulic fractures hf condc . The operational parameter considered in the shale gas development process is the bottom hole flowing pressure p wf . Augmenting all these parameters together, the optimization parameter set is: To initialize the ensembles of the optimization parameter set, 50 initial ensemble realizations are generated by sampling the parameters in the synthetic shale gas production model with a uniformly distributed random perturbation.
The parameter values used in the synthetic model and the optimization process are listed in Table 4.
During the optimization process, the EDFM-based shale gas production model is used to compute the 10-year cumulative shale gas production. The computation is conducted for each ensemble member in the generated initial parameter set, and the NPV is evaluated with considering the shale gas production revenue and the cost of the well drilling and completion. The covariance between the prior optimization parameters themselves, xx , and the covariance between the prior optimization parameters and the NPV values, xg , are computed by Eqs. (11) and (14), respectively. Then the optimization parameters in the initial ensembles can be updated by Eq. (16). The stopping criterion is that the increment rate of the NPV in two consecutive iterations is < 1%. Figure 6 depicts the optimization process of the design parameter for fractured horizontal well in the shale gas reservoir with natural fractures. Here, the iteration step is 33, which means that the optimization process stops after 33 (28) , , D, L, hf stage , hf xf , hf condc , p wf iterations. It can be seen that the design of the fractured horizontal well has been changed dramatically. Figure 6a visualizes the initial design of fractured horizontal well with the design parameters shown in the Table 4. Under the initial condition, the associated cumulative shale gas production is 4.57 × 10 7 m 3 and the NPV value is 34.84 million CNY. Figure 6b shows the optimization results after 4 iterations. The parameters ( , , D) characterizing the horizontal well location have been changed to 315.16 m, 89.01°, and 562.69 m. Due to the existence of the natural fractures, the distribution of the natural distribution makes the shale reservoir anisotropic and the adjustment of the well location is to improve the well productivity. The length of the horizontal well, L, has been increased to 676.45 m. The number of the hydraulic fracture stages has been slightly decreased to 10, but the half-length of the fracture has been increased to 184 m and the conductivity of the hydraulic fractures has been increased to 427.42 mD m. These parameter changes show that the optimization algorithm tends to find the most economical scenario to design the fractured horizontal well and improve the productivity. The bottom hole flowing pressure value has been decreased to 25.7 MPa. After four optimization iterations, the corresponding cumulative shale gas production is 6.46 × 10 7 m 3 and the NPV value is 61.52 million CNY. Figure   After 25 optimization iterations, the corresponding cumulative shale gas production is 9.42 × 10 7 m 3 and the NPV value is 115.87 million CNY. Figure 6f shows the optimization results after 33 iterations. The parameters ( , , D) characterizing the horizontal well location have been changed to 354.66 m, 78.36º, and 528.15 m. The length of the horizontal well, and the stage number, half-length and the conductivity of the hydraulic fractures have been increased to 817.2 m, 17, 229 m, and 621.4 mD m. The bottom hole flowing pressure has been decreased to 19.02 MPa. After 33 optimization iterations, the corresponding cumulative shale gas production is 10.01 × 10 7 m 3 and the NPV value is 143.42 million CNY. During the entire optimization process, it can be seen that the length of the horizontal well, and the stage number, half-length and the conductivity of the hydraulic fractures have been increased, and the bottom hole flowing pressure has been decreased. All the adjustments are conducive to the improvement of well productivity. In the meanwhile, the associated location adjustment of fractured horizontal well tends to maximize the connections with the natural fractures so that they can form an effective fracture network. It creates an advantageous condition for the shale gas to flow from the matrix to the wellbore.
The objective of the optimization is to maximize the NPV value through the obtained 10-year shale gas cumulative production from the simulation model with a given fractured horizontal design. Figure 7 shows the changes of the cumulative production of shale gas and the NPV during the entire optimization process. Both the shale gas production The cumulative production changes from its initial value 4.57 × 10 7 m 3 to 10.01 × 10 7 m 3 after optimization, and the NPV value changes from its initial value 34.84 million CNY to 143.42 million CNY after optimization. Since the objective function is NPV, it can be found that the NPV has been increased to 4.1 times more than its original value after optimization.

Conclusions
An ensemble-based optimization method is proposed to cope with the optimization design problem of the hydraulically fractured horizontal well in the shale gas reservoir. When building the simulation model of shale gas reservoir, the effects of gas-water two phase flow, gas absorption/desorption, and Knudsen diffusion in the shale matrix are considered. The equivalent permeability tensor method is used to consider the effect of the large number of developed natural fractures, which greatly improves the calculation efficiency. The embedded discrete fracture (EDFM) is used to explicitly account for the geometry of the hydraulic fractures so that the fracture flow can be better characterized. Based on the steepest ascent method, the ensemble optimization method is used to approximate the update gradient by covariance, which greatly reduces the difficulty of obtaining the update gradient. The well location, well length, number of fracturing stages, fracture half length, fracture conductivity, and bottom hole flowing pressure of fractured horizontal wells in shale gas reservoir are optimized simultaneously to achieve a global optimal performance. The horizontal well is parameterized through Hough transform method, which improves the continuity of the parameters and the Gaussianity required by the ensemble optimization algorithm. The optimization results show that the ensemble optimization algorithm is effective, and the NPV can be greatly increased to approach to its maximal value after optimization.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.