Horizontal drilling and multistage-hydraulic fracturing have been widely applied in the development of unconventional reservoir oil and gas resources. It is well known that many factors determine the productivity of horizontal wells, including formation petrophysical properties, fluid properties, and completion and fracturing parameters^[1]. Among them, formation petrophysical properties and fluid properties are uncontrollable factors, and completion and fracturing parameters are controllable. From the perspective of reservoir engineering, the aim of horizontal well fracturing is to increase contact area between fractures and reservoir, conductivity in propped hydraulic fractures, and reduce flow resistance to enhance well productivity. In the design of development well-pad with multiple fracturing wells, the main indexes need to be considered are well spacing, fracturing stage spacing, amount of proppant per fracture and fracture dimensions (length, width and height)^[2,3,4].

The difficulty in optimization of well-pad design is to make sure that the simulation can reflect the characteristics of pro-duction performance under the effects of inter-fracture interference, interwell interference, and fracture conductivity, so the relationships between main factors and production performance and the correlation between the main factors can be made clear. In the terms of performance simulation of multiple fracturing horizontal well (MFHW), several analytical and numerical methods have been developed. On the basis of instantaneous functions presented by Chen and Raghavan^[6], Wang et al.^[5] obtained the transient response of single MFHW by integrating the analytical solution for uniform-flux fracture with an impact function of fracture conductivity. Chen et al.^[7] proposed a semi-analytical method derived from an analytical reservoir solution and a numerical fracture solution to forecast the full life-cycle performance of MFHW, which has the flexibility of capturing the geometry complexity of fracture network and the interplay between reservoir and fractures. Fang et al.^[8] established a numerical model considering the multi-scale flow in discrete fractures, and simulated the production performance of a well pad with three volume frac-tured horizontal wells with long lateral section. Considering the well interference through both matrix permeability and hydraulic fracture hits, Yu et al.^[9] developed a numerical, compositional model in combination with an embedded discrete fracture model (EDFM) to simulate both the response to interwell interference, to find out the impact of well spacing on pressure interference in the adjacent wells.

In the aspect of analyzing the impact of main control factors, single factor control method is commonly used to examine the factors one by one, which ignores the correlations between the factors and thus could lead to wrong results^{[10,11,12,13]}. Although orthogonal test and grey correlation methods improve in limitedness than the single factor control analysis method, they are still multiple-scheme comparison methods, unable to cover all possible optimal solutions^[14,15,16]. Automated optimization algorithms such as artificial neural network and genetic algorithm can provide globally optimal design, but huge in computation, they are time-consuming and can’t be used widely, especially when the number of wells and fracturing stages are large^[17,18].

Systematic optimization of completion parameters of MFHWs is a complicated nonlinear problem in the field of reservoir engineering and mathematical optimization. On the basis of a semi-analytical modeling, this paper developed a comprehensive optimization workflow to maximize NPV by integrating extended proppant index method and globally optimal algorithm. The principle objective of this work is to provide a quantitative assessment of optimal decisions such as number of wells, number of fractures per well, mass of proppant and fracture dimensions.

The configuration of the MHFW system used in this study is shown in Fig. 1. All MFHWs are assumed to be evenly located and completed in an isotropic, homogeneous and horizontal formation in rectangular shape. The formation fully penetrated by hydraulic fractures has uniform thickness,permeability and porosity, all of which are constant. All transverse fractures are of identical properties. Besides, all MFHWs are produced at constant bottomhole pressure (BHP).

The fluid flow in the reservoir with low permeability is governed by non-Darcy’s law. Based on Ertekin model^[19], the flow velocity in the reservoir is defined as the sum of the component of the Darcy flow velocity and the slip flow velocity:

where the apparent gas slippage term is not a constant but a variable given by ${{b}_{\text{am}}}\text{=}\frac{{{D}_{\text{g}}}{{\mu }_{\text{g}}}{{c}_{\text{g}}}{{p}_{\text{m}}}}{{{K}_{\text{m}}}}\frac{1}{{{p}_{\text{m}}}}$.

The slip velocity is presented by modifying Fick’s law, which is caused by the concentration gradient in the matrix. At matrix-fracture interface, the concentration gradient is negligible^[20], and the slippage effect approaches zero, that is:

To linearize the gas-flow governing equations, the pseudo- function approach is defined. Pseudo-pressures in the fracture and reservoir are respectively defined by

And the nonlinear factor, γ_m, is introduced to account for the influence of gas slippage, which is expressed as

Correspondingly, the concept of pseudo-time in the reservoir is defined as

where β_m is a dimensionless rescaling factor which describes the behavior the fluid viscosity and compressibility during depletion, and λ_m is a dimensionless viscosity-compressibility ratio, which are defined as follows:

${{\beta }_{\text{m}}}(t)=\frac{1}{t}\int_{0}^{t}{{{\lambda }_{\text{m}}}(\tau )\text{d}\tau }$, ${{\lambda }_{\text{m}}}(t)=\frac{{{\gamma }_{\text{m}}}(t)}{\frac{{{\mu }_{\text{g}}}(t)}{{{\mu }_{\text{gi}}}}\frac{{{c}_{\text{t}}}(t)}{{{c}_{\text{ti}}}}}$

Substituting pseudo-pressure and pseudo-time into the gas governing equation renders gas diffusivity equation amenable to the liquid diffusivity equation.

Since the fluid flow caused by different MFHWs is confined to an identical pressure system, the fluid flow can be decomposed into continuous flow in fractures and flow in reservoirs. In the modeling of the fluid flow, two sets of separate coordinates were adopted for the fracture and reservoir, and the pressure and influx distribution along fracture face can be achieved by coupling the reservoir flow model and fracture model on the basis of continuity condition that the reservoir pressure equals to the fracture pressure on the fracture-reservoir interface.

To reflect the effect of interwell fracturing interference on the opening and propagation of hydraulic fractures, a fracture model with variable width is adopted in this study. To our knowledge, when the propagation direction of hydraulic fracture is parallel to the maximum stress direction, the farther from the wellbore, the stronger the interwell stress will be, thus the net pressure to open up fracture will increase, the migration of proppant will be blocked, and the fracture width will decrease gradually from the wellbore outward^[21]. That is expressed as:

For simplicity, dimensionless pseudo-pressure in reservoir, dimensionless pseudo-pressure in fracture, dimensionless production rate and dimensionless cumulative production are respectively defined by

and dimensionless time, dimensionless length, dimensionless conductivity and dimensionless influx density are defined as follows:

where χ and ψ are scaling factors. In SI unit, χ=0.5/π and ψ=1; in Darcy unit, χ=1.842 and ψ=0.0036.

Since all MFHWs are located in a pressure system, the drainage area can be divided into a set of sub-drainage areas according to the principle of pressure superposition. Fig. 2 shows the equivalent sub-drainage area drained by single fracture (sub-drainage area equals to 1/n_wn_f of total drainage area). Using pseudo-function definitions, the equation governing gas flow in the reservoir is given by

where S_f(x,y,t) is the source function for influx distribution of hydraulic fracture with the length of L_f:

To account for the effect of gas desorption, q_de(t) is defined to describe the supply of gas desorption^[7]:

The hydraulic fracture can be regarded as the Green’s function for a line source along which flux density is satisfied as q_fD(s). The gas desorption is considered in the dual-porosity model of f(s). After using the Laplace transformation, Equation (16) can be written as the following dimensionless form,

From the perspective of reservoir flow model, the fracture can be discretized into a set of segments, so Eq. (19) is rewritten into a discretized equation in which each segment has a length $\Delta {{x}_{fDi}}$and flux q_fDi. Based on the principle of pressure superposition, the pressure at any point in the fracture can be expressed as:

where $\tilde{p}$_uDj,i represents the pressure disturbance to the j-th segment caused by i-th segment. According to the Green’s function^[22], the dimensionless pressure solution for uniform- flux segment can be written as

where

${{\tilde{H}}_{n}}=\frac{\cosh \left[ {{\varepsilon }_{n}}({{y}_{\text{eD}}}-|{{y}_{\text{D}j}}-{{y}_{\text{wD}i}}|) \right]}{\sinh ({{\varepsilon }_{n}}{{y}_{\text{eD}}})}\text{+}\frac{\cosh \left[ {{\varepsilon }_{n}}({{y}_{\text{eD}}}-|{{y}_{\text{D}j}}+{{y}_{\text{wD}i}}|) \right]}{\sinh ({{\varepsilon }_{n}}{{y}_{\text{eD}}})}$${{\tilde{F}}_{ji}}=\int_{-\text{ }\!\!\Delta\!\!\text{ }{{x}_{\text{fD}i}}}^{\text{ }\!\!\Delta\!\!\text{ }{{x}_{\text{fD}i}}}{\left( \frac{2\pi }{{{x}_{\text{eD}}}}\sum\limits_{n=1}^{\infty }{\frac{1}{{{\varepsilon }_{n}}}\cos \frac{n\pi {{x}_{\text{D}j}}}{{{x}_{\text{eD}}}}\cos \frac{n\pi {{u}_{\text{D}}}}{{{x}_{\text{eD}}}}} \right)d{{u}_{\text{D}}}}$${{\beta }_{\text{wD}j}}=({{x}_{\text{D}j}},{{y}_{\text{D}j}})$ ${{\varepsilon }_{n}}=\sqrt{sf(s)+{{n}^{2}}{{\pi }^{2}}/x_{\text{eD}}^{2}}$

$f(s)=\omega +\frac{{{a}_{\text{g}}}{{D}_{\text{D}}}(1-\omega )}{s+{{D}_{\text{D}}}}$

For simplicity, Equation (20) is rewritten as the corresponding matrix:

In 1D coordinate of fracture model, the fluid flow from the reservoir into hydraulic fracture, and then flow towards wellbore along fracture panel. As shown in Fig. 3a, the flux density along fracture is given by q_f(x,t), and the rate extracted from wellbore is denoted as q_w(t).

To account for variable conductivity fracture shown in Fig. 3b, the pseudo-pressure equation is described in the following dimensionless form:

where δ is the Dirac function, and dimensionless variable conductivity is written as

Equation (23) is a pseudo-linear differential equation because the conductivity C_fD is a function of x_Dhf. Here coordinate transformation equation is introduced to make the pseudo- linear problem linear, which is expressed as

and ${{\hat{C}}_{\text{fD}}}=\frac{{{L}_{\text{fD}}}}{\int_{0}^{{{L}_{\text{fD}}}}{C_{\text{fD}}^{-1}\left[ {{p}_{\text{fD}}}({{x}_{\text{D}}}) \right]d{{x}_{\text{D}}}}}$. Incorporating Equation (25) and boundary element method^[23] to deal with Equation (23), the solution for variable conductivity fracture can be derived as

where G is the integral of Heaviside unit step function, and the Fredholm integral function is defined as

To take the additional pressure drop due to the convergence of fluid into the horizontal wellbore from transverse fractures, the convergence flow skin S_c is introduced to modify Equation (26)^[24], which is rewritten as

${{S}_{\text{c}}}=\frac{2{{h}_{\text{D}}}}{{{{\hat{C}}}_{\text{fD}}}{{L}_{\text{fD}}}}\left[ \ln \left( \frac{{{h}_{\text{D}}}}{2{{r}_{\text{wD}}}} \right)-\frac{\pi }{2} \right]$

Equation (28) can be also expressed in the matrix form:

where ${{\gamma }_{i}}=G({{\xi }_{Di}},{{\xi }_{wD}})$, ${{\alpha }_{j,i}}=R{{S}_{i}}+R{{T}_{i}}$,

$G(x,{{x}_{0}})=\left\{ \begin{align} & x-{{x}_{0}}\ \ (x\ge {{x}_{0}}) \\ & 0\ \,\ (x<{{x}_{0}}) \\ \end{align} \right.$,

$R{{S}_{i}}=\left\{ \begin{align} & {{\xi }_{\text{D}j}}\text{ }\!\!\Delta\!\!\text{ }{{\xi }_{\text{D}j}}-0.5(\xi _{\text{oD,}i}^{2}-\xi _{\text{oD,}i-1}^{2})\quad (i<j) \\ & {{\xi }_{\text{D}j}}({{\xi }_{\text{D}j}}-{{\xi }_{\text{oD,}j-1}})-0.5(\xi _{\text{D}j}^{2}-\xi _{\text{oD,}j-1}^{2})\quad \ (i=j) \\ \end{align} \right.$,

$R{{T}_{i}}=\left\{ \begin{align} & ({{\xi }_{\text{wD}}}-{{\xi }_{\text{D}i}})\text{ }\!\!\Delta\!\!\text{ }{{\xi }_{\text{D}j}}\quad (i<1+\Xi \text{)} \\ & \left( {{\xi }_{\text{wD}}}-\sum\limits_{i=1}^{\Xi }{\text{ }\!\!\Delta\!\!\text{ }{{\xi }_{\text{D}i}}} \right)\left( 0.5{{\xi }_{\text{wD}}}+\sum\limits_{i=1}^{\Xi }{\text{ }\!\!\Delta\!\!\text{ }{{\xi }_{\text{D}i}}} \right)\quad (i=1+\Xi ) \\ \end{align} \right.$

According to the continuity condition that the reservoir pressure equals to the fracture pressure on the fracture-reservoir interface, the relationship is expressed as:

It is worth noting that the relation between flux-density of fracture segment and rate of wellbore in the transformed domain is:

where${{f}_{\text{c}j}}=\Delta {{\xi }_{\text{D}j}}/\Delta {{x}_{\text{D}j}}$. The element of $\tilde{p}_{\text{uD}j,i}^{{}}({{\beta }_{\text{wD}j}},{{\beta }_{\text{wD}i}})$ in the matrix of Equation (22) is modified as${{f}_{\text{c}j}}\tilde{p}_{\text{uD}j,i}^{{}}({{\beta }_{\text{wD}j}},{{\beta }_{\text{wD}i}})$.

According to coupled condition [Equation (30) and Equation (31)], combining the solution for the reservoir [Equation (22)] and the solution for the fracture [Equation (29)] provides the closed solution for flux-density distribution along fracture and production rate of the well. The unknown variables can be expressed in vector form:

The cumulative production is obtained by using Laplace transformation, which is given by $\tilde{G}$_pD=$\tilde{q}$_wD/s.

It is worth noting that the pseudo-time is a function of spatial and temporal variables. The pseudo-time can be regarded as single integral of average reservoir pressure over time, as presented in Equation (33), and can be calculated by incorporating material balance equation [Equation (34)]. The detailed information on iteration algorithm is provided in Reference [25].

Afterwards, to verify the accuracy of the semi-analytical model in our work, the transient response solutions were compared with other semi-analytical model and the results of numerical simulation. For the semi-analytical results presented by Poe et al.^[26], the variation of fracture conductivity is described by four constant-length zones and each zone is assigned a different conductivity. As shown in Fig. 4, our pressure solutions are consistent with Poe’s results, validating the calculation accuracy of our model in simulating pressure in one fracture with variable width.

Commercial numerical simulator was used to further verify the reliability of our model in simulating interwell interference. The parameters used in the model were: p_i=45.0 MPa, p_wf=6.5 MPa, h=20 m, s_o=85%, c_t=4.35×10^-4 MPa^-1, K_m= 0.01×10^-3 μm²; μ_o=1 mPa·s, B_g=1.0 m³/m³, r_w=0.019 m; n_f=30, x×y×z=1 500 m×1 600 m×20 m, x_f=114.95 m, w_f=0.0127 m, ϕ_m=35%. To generate accurate rate response, the local-grid-refinement (LGR) was used to refine the grids around fracture panel. Since φ_f=35% and w_f=0.0127 m, the equivalent fracture porosity calculated is 4.4%, the equivalent fracture width is 0.1 m, and the equivalent fracture permeability is calculated using equation K_fmax=12C_fd, which equals to 37.68×10^-3, 376.8×10^-3, 3 768×10^-3 respectively in this study, the permeability distribution is presented in Fig. 5.

Due to the ideal assumption of even well configuration, all the wells have the same pressure distribution. Fig. 6 shows the dimensionless production performance from numerical simulation and the method presented in this paper, it can be seen the single well production and cumulative production from the two methods tally well with each other, verifying the adaptability of the method proposed in this paper in cases with multiple wells and multiple fractures.

As mentioned before, the production performance of well pad is affected by various factors, including fracture dimensions (width, length and height), number of wells, number of fractures, and fracture/well arrangement. To simplify the issue, we assume an ideal homogeneous configuration, with wells and fractures in even distribution and same in dimensions. The fractures occupy the same sub-drainage area (1/(n_wn_f) of the total drainage area). Thus the factors affecting production performance of the well pad can be reduced to length and conductivity of the fracture.

At a given mass of proppant, the fracture length and fracture width are competing for the proppant. The maximum productivity index would be achieved on the best compromise between length and conductivity. The classical proppant index method theory is originally proposed for single fracture with uniform conductivity in the pseudo-steady state (PSS) flow conditions^[27,28], and the objective function is the productivity index in the PSS condition, while not (cumulative) rate in transient conditions.

In this paper, an extension of proppant index method theory to account for the geometry of multiple MFHWs is presented, and the proppant number is modified as

The extended proppant index method could consider time dependency of transient response covering both transient and PSS conditions. The fixed proppant number in the proppant index method is regarded as technical constraint. With further consideration of economic constraint, a nested optimization method containing outer-optimization and inner-optimization shells is further developed to optimize an overall economic objective.

Fracture length and conductivity (set the fracture height is the same as the formation thickness) are key parameters of inflow and outflow relations between fracture and reservoir. Considering about the positive relationship of fracture width and fracture permeability, fracture conductivity ability (K_fw_f) is used to represent fracture width. We take a single fracture as an example to discuss the optimization process of fracture length and conductivity. The fracture size are assumed as 400 m×50 m×20 m, where the fracture length is represented by fracture-penetration ratio (I_x=L_f/x_e), and the fracture width is represented by dimensionless conductivity.

The optimization results of fracture conductivity and length with proppant index in Fig. 7 show when proppant volume is not considered, estimated ultimate recovery (EUR) denoted by dashed line increases monotonously with the increase of I_x and C_fD till reaches the maximum value. At this point, the fracture full penetrates the drainage area and reaches infinite conductivity. When the constraint of proppant volume is considered, the cumulative production and fracture dimensions have optimal values, which are denoted by discrete points in the Fig. 7. The maximum conductivity and length cumulative production increases with the increase of proppant index. More interestingly, the optimum fracture conductivity and length are different at different times as seen from Fig. 7a, 7b.

To determine the effect of production time on the optimum results of fracture conductivity and length, the optimized fracture dimensions with time were calculated (Fig. 8). As shown in Fig. 8a, as expected, EUR_max is an increasing function of t_D for N_prop=0.1, while C_fDopt is a decreasing function and approaches a limit. From Fig. 8b to Fig. 8d, we can see that the larger the proppant index, the larger the dimensionless cumulative production at the same time will be, the wider the variation range of optimum dimensionless conductivity in the same time interval, and the larger the value the dimensionless conductivity approaches will be. For example, as presented in Fig. 8b for N_prop=1, when t_D is in the range of 0.01-1000, the corresponding C_fD,opt is in the range of 1.62-100.00, and the limit value of C_fD,lim is approaching 1.62. In contrast, in Fig. 8d for N_prop=1000, the variation range of C_fD,opt is 15-400 with C_fD,lim=15.

Fig. 9 shows the variation patterns of maximum EUR and optimal fracture conductivity and length at different proppant indexes and different times. Fig. 9a shows a positive correlation between C_fD,opt and N_prop, which means that the C_fD,opt decreases with the decrease of N_prop till reaches a limited value and then keeps constant basically. The bigger the value of t_D is, the larger the variation range of N_prop. For example, when t_D=1000, C_fD,opt approximates a constant in the range of N_prop=1×10^-8-1×10^-1. However, when t_D=0.01, N_prop is in the range of 1×10^-8-1×10^-7. Besides, both C_fD,opt and t I_x,opt increases with I_x. Once I_x is fixed at unity, C_fD,opt would show a linear relation with N_prop. Fig. 9b shows that for a given N_prop, C_fD,opt decreases with time and approaches a constant value. The constant value of C_fD,lim is positively correlated with N_prop [e.g., C_fD,opt=1.62 for N_prop=1×10^-2, while C_fD,opt=175 for N_prop=1×10³]. Put it in another way, when dimensionless fixed period of time is shorter, the optimum design indicates a fracture with shorter length and higher conductivity; when dimensionless fixed period of time is longer, the opposite is the case. It is worth noting that the characteristic value of C_fD,opt is 1.62 in the case of small N_prop, which is consistent with the results of Valko and Economides^[27], proving the method presented in this paper is reliable.

Assuming the drainage area is fixed, the geometrical shape affects the well performance and the optimum fracture conductivity and length. The drainage-area aspect ratio describes the geometrical shape, which is defined as λ=y_e/x_e. Noted that aspect ratio is also determined by n_w and n_f (y_e=y_ef/n_f, x_e=x_ef/n_w). Fig. 10 shows the effect of λ on the production performance. In this case, the fractures are assumed to fully penetrate the pay in the plane, and all fractures have infinite conductivity. The maximum dimensionless accumulative production (0.5x_eDy_eD/π) could be reached when the production decreased to 0 with a fixed drainage area, irrelevant to the ratio of λ. The initial rate is higher and the maximum EUR is reached in a shorter duration for drainage area with more pronounced elongation when λ is smaller. In fact, the more elongated the drainage area, the more efficiently it can be drained, since the corresponding fracture length is larger, thus making it easy to drain the whole drainage area.

When the proppant number is taken into account, Fig. 11 shows the effect of λ on maximum EUR and the optimal fracture conductivity and length when t_D is small. Fig. 11a shows that EUR_max increases with N_prop until reaches the peak. The smaller the λ (small fracture space and large well space), the larger the maximum dimensionless cumulative production will be, and at this point, the fractures penetrate the whole formation and reach infinite conductivity (Fig. 11b). When the proppant index is smaller (less than 1×10³), for different λ values, the corresponding optimum conductivities are the same, but the optimum fracture-penetration ratio is larger with larger value of λ (Fig. 11b), so the corresponding maximum dimensionless cumulative production is larger (Fig. 11a). It must be highlighted that the production characteristics reflected by the template in Fig. 11 is related to the production time; when the time is long enough, the peak values of EUR_max for all λ values are equal; when the time is short enough, λ has little effect on type curves (shown in Fig. 11) in the case of smaller N_prop.

The development effect of multi-stage fracturing horizontal well is optimized by increasing the contact area between fractures and formation, reducing interwell and inter-fracture interferences, balancing the inflow and outflow between fracture and formation. When these seepage relationships reach balance, the production effect is the best. In this study, half of the well pad is assumed to be 1 500 m×1 500 m×20 m, with the proppant volume as constrain, and the NPV as the objective function, a nested method was used to optimize the parameters. The detailed workflow entails the following steps.

Step 1: Define input variables, including reservoir/well properties, fluid properties, proppant type and the targeted production period.

Step 2: Define three key variables to be optimized, including number of wells n_w, number of fracturing stages per well n_f, and proppant volume per fracture, M_p.

Step 3: Calculate the optimum fracture conductivity and length and maximum cumulative production at different number of wells, fracturing stages and proppant volumes by using the proppant index method.

Step 4: Calculate the corresponding NPV.

Step 5: Based on Powell global-optimization algorithm, return to Step 2 till the maximum NPV is achieved, and the horizontal well and fracturing parameters at this point are the best design parameters.

The NPV is calculated as

In this case, the selected dimensionless fixed period of time is set to be 1000, and total mass of proppant of multiwell pad is assumed to be fixed, as the actual fracturing scale is limited by engineering conditions. Fig. 12 shows a 3D plot of the total EUR_max variation with the combination of n_w and n_f. It can be seen from Fig. 12, EUR_max increases with the increase of wells and fracturing stages monotonously. When the proppant index is smaller, EUR_max gradually increases in the same magnitude with the increase of wells and fracturing stages (Fig. 12a). When proppant index is bigger, the increase magnitude of EURmax decreases relativley. When the fracturing stages are more than 30 and well pads are more than 5, the EURmax hardly increase anymore (Fig. 12b). Therefore, in the case with the EURmax as the objective function and the fracturing scale as the constraint, there is only optimal fracture conductivity and length but no optimal well space (number of wells) and fracture space (number of fractures) for a well pad.

Substituting the development index results into the economic evaluation model, the optimization was repeated with NPV as the objective function, the results are shown in Fig. 13. There are obvious extreme points in Fig. 13, indicating there are optimum well space and fracture space in this case. Similarly, the corresponding optimization results of fracture conductivity and length are shown in Fig. 14. It can be seen with the increase of n_w and n_f, although the EUR_max increases, the costs also go up. In other words, when the incremental NPV caused by EUR increase is offset by the cost increase caused by creating more wells and fractures, the economic benefit would get worse gradually. Therefore, under the dual constraints of fracturing scale and economic benefit, there are optimum n_w, n_f and fracture conductivity and length for a well pad, providing optimization room in development program design.

It can be seen from Figs. 13 and 14, when the proppant index is smaller (N_prop=0.1), the optimal well space is smaller (n_w,opt≈8), the optimal fracture space is wider (n_f,opt≈35), the corresponding optimal fracture-penetration ratio is smaller (I_x,opt≈11%), and the fracture conductivity is smaller (C_fD,opt≈ 1.60). When the fracturing scale is larger (N_prop=1×10⁴), compared with the case with small fracturing scale, the optimal well space increases (n_w,opt≈6), the optimal fracture space decreases (n_f,opt≈50), the fracture penetration rate increases significantly (I_x,opt≈72%), but the optimal conductivity is still fairly low (C_fD,opt≈3.54).

The specific results of nested optimization are related to input parameters, such as geometry of the well pad, development duration, economic parameters, complexity of fracture networks and constraint conditions (i.e., mass of proppant per fracture, well or well pad etc.). However, these general statements can serve as reference in the more-detailed optimization works.

In this study, we established a general semi-analytical model to study the production performance of a well pad with multiple MFHWs. Dimension transformation is incorporated to change the flow equation of the variable conductivity fracture into the equation of the constant conductivity fracture, and a universal approach is developed to solve the coupled reservoir-fracture flow model.

In the inner-optimization shell, there exists a best compromise between fracture conductivity and fracture length under the constraint of fracturing scale. As proppant index and time increase, the maximum accumulative production increases monotonously, the optimal fracture penetration rate also increases. With the decrease of proppant index and increase of time, the optimal conductivity of fracture decreases gradually and approaches to a limited value (C_fD,lim=1.62). Besides, the ratio of fracture space to well space also affects the results of the optimal fracture conductivity and length. In the case with larger N_prop, the design with larger well space and smaller fracture space is better.

When economic benefit isn’t considered, the EUR_max increases with the increase of n_w and n_f, there is no optimal wells and fracturing stages. But the EUR increase magnitude decreases with the increase of proppant index, when the n_w and n_f reach a critical value, the EUR doesn’t increase anymore.

In the nested optimization accounting for technical and economic constraints, there exists room for optimization of fracture design when the NPV increment can be offset by increase in cost. When the N_prop is relatively lower, the design with more wells and short-length fractures is suggested; when the N_prop is higher, the design with less wells and longer fractures is better.

a_g—gas desorption coefficient, dimensionless;

B_g—gas volume factor, m³/m³;

b_am— slippage factor, dimensionless;

C_fD—dimensionless fracture conductivity;

C_fracture— cost of fracturing one stage, Yuan;

C_well—the cost of drilling one well, Yuan;

c_t—formation compressibility factor, Pa^-1;

c_g—gas compressibility factor, Pa^-1;

D_D—dimensionless diffusivity coefficient;

D_g—gas diffusivity coefficient in matrix, m²/s;

FC—the total fixed cost, Yuan;

G_p—cumulative production, m³;

G_p,ii—the ii-th year cumulative production, m³;

h—formation thickness, m;

I_x—fracture penetration ratio;

ii—year number;

i_r—annual interest rate, dimensionless;

K_f—fracture permeability, m²;

K_m—reservoir permeability, m²;

L_f—fracture length, m;

L_ref—reference length, m;

m_f—pseudo-pressure in fracture, Pa;

m_m—pseudo-pressure in matrix, Pa;

N—fracture micro-stage;

N_prop—proppant index;

NPV—net present value, yuan;

n—infinite series;

n_f—number of fractures per well;

n_w—number of wells in a well pad;

n_year— production time, a;

OGIP—original gas in place, m³;

P_gas—gas price, Yuan/m³;

p—pressure, Pa;

p_f—initial pressure of fracture, Pa;

p_m—formation pressure, Pa;

p_w—wellbore pressure, Pa;

p_wf—bottomhole pressure, Pa;

q_f—influx density, m²/s;

q_m—formation flow density, m³/s;

q_w—production rate, m³/s;

$\tilde{q}_{\text{fD}n,i}^{\quad \left\langle w \right\rangle }$—dimensionless production of the i-th segment of n-th fracture in w-th well;

R_m—pore radius of gas desorption, m;

r_w—wellbore radius, m;

s—Laplace variable, dimensionless;

Sc—conflux skin factor, dimensionless;

T—Temperature, K;

t—time, s;

t_a—pseudo-time, s;

V_prop—volume of proppant, m³;

V—gas concentration, m³/m³;

V_E—gas concentration in the center of spherical pore, m³/m³;

v_m—flow velocity in matrix, m/s;

w_f—fracture width, m;

x, y—formation system coordinate, m;

x_e—width of sub-drainage area, m;

x_ef—reservoir width, m;

x_f—fracture half-length, m;

Δx_f—fracture micro-stage length, m;

x_hf—coordinate in fracture, m;

x_o, y_o—end points of fracture, m;

x_whf—intersection between fracture and wellbore, m;

y_e—length of sub-drainage area, m;

y_ef—reservoir length, m;

Z_g—gas deviation factor, dimensionless;

Z—deviation factor, dimensionless;

β_m—dimensionless rescaling factor;

γ_m—nonlinear modified factor;

δ—Dirac function;

ζ, ξ, ς, u, τ—general spatial variable;

λ—ratio of fracture spacing and well spacing;

λ_m—dimensionless viscosity-compressibility ratio;

μ_g—gas viscosity, Pa·s;

ξ_D—dimensionless transformed spatial variable;

ϕ_m—reservoir porosity, %;

ψ—unit conversion factor;

ω—storage ratio, dimensionless;

—micro-element number of fracture;

▽p_m—pressure gradient, Pa/m.

avg—average;

D—dimensionless;

i—initial pressure condition;

i, j—number of fracture segment;

opt—the optimal;

max—the maximum;

sc—standard condition;

w—wellbore position.

~—Laplace symbol.