A new approach to finding effective parameters controlling the performance of multi-stage fractured horizontal wells in low-permeability heavy-oil reservoirs using RSM technique

The application of multi-stage fractured horizontal well (MSFHW) due to its costly operation necessitates optimization of associated fracture parameters to ensure its economic success. In comparison to significant number of studies dedicated to use of MSFHWs for shale gas reservoirs, there are only few researches available for oil systems. This study explores the optimum criteria for a number of important fracture parameters in low-permeability heavy-oil systems. For this purpose, a response surface methodology (RSM) was employed to examine the simultaneous effect of four fracture parameters, including the number of fracture stages, fracture length, fracture width and fracture conductivity, on well productivity. The evaluations were conducted on two homogeneous and heterogeneous permeability systems. The optimization of fracture parameters was also performed on an economic basis by utilizing the net present value (NPV) concept. Useful charts were also generated providing practical insights into the individual and combinational effects of fracture parameters on well performance. The results from this study demonstrated that the fracture conductivity and the number of fracture stages were, respectively, the first two important parameters controlling the well productivity for rock systems with higher permeability. However, when rock texture became tighter, the number, and to a lesser extent the length, of fractures exhibited more evident role on production improvement, especially in the case of heterogeneous reservoirs. The results also underlined the significance of economic considerations, in particular, when determining the optimum fracture length and number of fracture stages.


Introduction
Carbonates, because of their low-permeability rock textures, are good candidates for acid fracturing jobs. In this method, acid under high pressures is injected into the target zone 1 3 during which the surface of the rock minerals is dissolved unevenly. As a result, a non-uniform pattern is etched on the rock surface, which keeps the fracture tip open constantly, even after the pressure is released. Horizontal well drilling and hydraulic fracturing are two recently used IOR methods utilized to improve the well productivity owing to increased well-to-reservoir contact. Obviously, there is no question left why combining these two techniques, i.e. multi-stage hydraulic fracturing of horizontal wells (MSFHW), has gained significant attention for development of hard-toproduce resources. However, as the MSFHW technology is a complex and costly process, optimization of the fracture parameters is a crucial task to ensure an economic operation and satisfactory well performance. There are various methods developed for modeling the horizontal wells with multiple transverse fractures (Deng et al., 2014;MoradiDowlatabad and Jamiolahmady 2018;Raghavan and Joshi 1993;Valko and Amini 2007;Wang et al. 2010;Zhao et al. 2016). These models can be used to predict the performance of a fractured well Wang et al. (2016) combined response surface methodology with a flow-stress-damage model to estimate the stimulated reservoir area as a function of several parameters such as dip angle, stress difference and injection rate. They concluded that injection rate is a positive factor increasing the stimulated area since at higher rates less leakage into the formation discontinuities happens. Tang et al. (2018) investigated the interactive effect of heterogeneity with some factors like in-situ stress gradient and stress shadow. Their results underlined the significance of considering formation heterogeneity in designing the fractures. Hareland and Rampersad (1994) developed a pseudothree-dimensional hydraulic fracturing model to optimize the controllable fracturing parameters such as pump rate, fluid rheology, proppant schedule and fracture length. Their results demonstrated that NPV does not necessarily increase with increasing the fracture length. Jabbari and Benson (2013) performed a case study for the optimization of hydraulic fracturing design in a low-permeability shale reservoir. Their results demonstrated the necessity of creating longer fractures in low-permeability rocks. Jahandideh and Jafarpour (2016) studied the impact of heterogeneity in shale brittleness and fracability on fracturing process and emphasized that the optimal number, location and length of fractures are dependant on rock fracability distribution. Romero et al. (2003) used a direct boundary element method to calculate the well productivity index of fractured wells. They found that the optimum range for fracture conductivity and its geometry is significantly affected by the fracture face skin. Besler et al. (2007) studied the fracture performance in horizontal wells and concluded that while longitudinal fractures exhibit low-conductivity requirements, increase of productivity in wells with transverse fractures is dependent on increase in fracture conductivity and width. Rahman et al. (2014) investigated the optimization of fracture configuration under geomechanical considerations. They concluded that elongating the fracture is more beneficial than increasing the fracture conductivity. Soliman and Grieser (2010) studied the effect of fracture spacing and timing of fracturing. They showed that the optimum fracture spacing is influenced by stress interference between fractures and fluid-flow conditions. Yang et al. (2017) made a study on optimization of the number of fracture stages and fracture spacing. They concluded that by earning a better description of reservoir heterogeneity, it is not anymore imperative to space fractures evenly; instead, they can be placed in intervals with higher density of natural fractures.
Most of the previous works on optimization of fracture design for multi-stage fractured wells focus mainly on single fracture parameters and their individual effect on well performance. This work presents a full analysis of the simultaneous effect of multiple fracture parameters governing the performance of MSFHWs and the level of importance of each parameter for various reservoir cases. A response surface methodology (RSM) has been used for this purpose. Fracture length, fracture width, fracture conductivity and number of fracture stages are four main parameters considered for two homogenous and heterogonous reservoirs. The net present value (NPV) has also been used to find the optimum fracture design criteria from an economic point of view. A description of the reservoir models, the range of fracture parameters and the methodology used for analyzing the results are described in the following sections.

Description of reservoir model
A single-well sector model of a carbonate reservoir located in the south-west of Iran was used for the main parts of this study. In addition, three other permeability models described later were imposed on the selected sector model to further investigate the sensitivity of the MSFHW performance to rock permeability. The oil bearing layer has a thickness of 250 ft, where a horizontal well, 3000 ft long, is penetrated through its middle. The no-flow-boundary reservoir has an initial pressure of about 4700 psi and contains a relatively heavy-oil fluid with an API of 22.4 and bubble point pressure of 1663 psi. The single well model produces under a constant bottom-hole pressure of 2000 psi.
The Eclipse commercial software was used to simulate the reservoir model (Eclipse Simulation Software 2010). An explicit design was implemented to model the fracture geometry by ascribing fracture properties including fracture length, width and permeability to the corresponding grid blocks around the wellbore. The fractures were considered to be transverse to the horizontal well and propagated across the full reservoir thickness. To accurately model the flow behavior around the fracture face, local grid refinement (LGR) technique was employed as shown by Figs. 1 and 2. The simulation model in Cartesian system consists of 177 × 35 × 5 grid blocks, respectively, in x, y and z directions. The final dimensions of the sector reservoir model are listed in Table 1.
As mentioned earlier, the main focus of this study is on the performance of MSFHW and level of importance of fracture parameters in a low-permeability heavy-oil reservoir (Case 3 in Table 2), which is known as the base case hereafter. The average horizontal permeability of this sector model is 4.7 mD. In addition, to complement this work and to examine the effect of fracture parameters with respect to rock permeability, three other permeability models were also created. These three cases, as listed in Table 2, include two homogenous models with permeabilities of 5 mD (Case 1)   and 0.5 mD (Case 2) and a heterogonous reservoir with an average permeability of 0.47 mD (Case 4). It should be noted that properties of these three cases, except their permeability and porosity, are similar to Case 3 (base case). The permeability map and its distribution in each layer of Case 3 are depicted by Figs. 3 and 4, respectively. More details of the properties of the reservoir model used in this study can be found in Table 3.

Methodology
Four fracture parameters including fracture length, fracture width, fracture conductivity and the number of fracturing stages were selected for this study. The range of these parameters will be discussed in the next sections. In total, 600 simulations were conducted for each reservoir model. The MINITAB software was then used to analyze the results of the numerical simulations performed (Minitab Software 2018). Accordingly, a response surface methodology (RSM) was employed to scrutinize the effect of fracture parameters on well productivity. An evaluation exercise was also carried out on the cost of fracturing operation to find the optimum fracture parameters from an economic perspective.

Selected range for fracture parameters
Al-Ameri and Gamadi (2019) in their study demonstrated that it is possible to have a full combination of fracture length, fracture width, and fracture conductivity by controlling the acid properties and the post-flush stage. Controlling the operational parameters such as pumping rate, injection rate, acid type, injection time and fluid loss control will help the engineers to adjust the fracture parameters and hence to maximize the fracture conductivity near the wellbore (Aljawad et al. 2019). According to studies performed on hydraulic fracturing and real practices carried out on fracturing operations in the oil and gas industry, the appropriate range for fracture parameters were selected for the purpose of this study as follows.

Fracture length
Emulsified acids are capable of creating longer fracture length without deterioration of fracture conductivity (Navarrete et al. 1998). Alghamadi (2006) expressed that the effectiveness of a hydraulic fracturing operation is strongly dependent on the effective fracture length and conductivity.
Under low temperatures, the fracture propagates to short distances due to inadequate stimulation. At high temperatures also the etched fracture length is limited because of the high reaction rate. All these parameters are controlled by various additives such as surfactants, retarders and emulsified acids (Dehghani et al. 2019). Williams and Nierode (1972) studied the acid penetration distance and stated that the fracture geometry is influenced by formation temperature, acid injection rate, acid concentration, and rock type. Creating a desired fracture length is achieved by designing the desired acid fluid system in many cases. Acid fluid loss during acid-fracturing treatments significantly limits the effective fracture lengths obtained in carbonate reservoirs. The fluid-loss control systems like a gelled-acid system (White et al. 1992) are effective on controlling the fracture length. Mukherjee and Cudney (1993) introduced a new gelbased acid that eliminates wormhole growth and increases the length of the fracture. The viscosity of the acid is temporarily increased from 30 to 1000 cP after the acid is injected into the formation, preventing acid leakage and leading to create more extended fracture length with the same volume of injected acid. Zhou et al. (2007) developed leak-off control acid (LCA) systems to reduce acid leak-off. The initial viscosity of LCA designed at a low value .The viscosity of LCA increased impressively when the value of pH decreases to a  specific value, and this can reduce the leak-off and enhance acid penetration. Rahman (2010) investigated the effect of injected acid rheology represented by the power-law exponent, injected rate and initial oil flow rate and interaction between them on fracture length and etched fracture width. He showed that fracture length decreases and etched fracture width increases by increasing power-law exponent.
de Rozieres et al. (1994) investigated the diffusion coefficient to predict the fracture length. They concluded appropriate design for acid fluid by a suitable combination of different acid type and diffusion coefficient, can significantly affect the fracture length and subsequently, the production of the fractured well.
Specific equipment can be used to provide the desired fracture length during acid fracturing operations. For example, diverters are commonly used along the stimulated wellbore to create more uniform stimulation and to enhance acid penetration radius (Aljawad et al. 2019). In summary, creating an appropriate fracture length is related to a number of factors, including the volume of injectable fluid, flow rate and fluid injection pressure, the type and viscosity of the injected fluid, and the amount of fluid that leaks into the formation during the fracturing process (Economides and Nolte 2000). Many of these factors can be controlled by choosing an appropriate design and selecting a suitable fluid system.
According to the real data collected from the fracturing operations in different fields, to cover a reasonable range, five fracture half-length (X f ) values from 50 to 400 ft was used in this work.

Fracture width
Fracture aperture or width (W f ) is one of the chief parameters affecting the fracture performance. For propped fracturing treatments, desired fracture width can easily be reached by appropriate proppant pack selection. However, in acid fracturing, initial fracture width is a function of the rate of acid injection, acid type, acid pumping time and acid volume, acid contact time, acid temperature and etching pattern (Penaloza 2013;Pournik 2008;Suleimenova 2015). It has also been shown that the fracture width increases with decreasing the number of fracture stages (Al-Ameri and Gamadi 2019). An appropriate fracture width can be achieved by controlling these factors.
For the reservoir model used in this study, four explicit values for fracture width ranging from 0.09 to 0.36 inches were selected.

Fracture conductivity
Fracture conductivity (C f ) is the product of fracture width (W f ) and fracture permeability (K f ) (Eq. 1). Suleimenova (2015), in his experimental work found out that not only the fracture permeability, but also the cross-section area for flow or fracture width, are crucial parameters affecting the flow capacity of the fracture. Acid fracture conductivity can be controlled by choosing proper operational parameters, although reaching a definite value is not assured. It has been demonstrated that extremely high contact times and high operation temperature do not necessarily guarantee higher fracture conductivity (Melendez Castillo 2007).
Acid properties such as viscosity, concentration and type affect the fracture conductivity. Also, the post-flush stage has an essential role on the magnitude of fracture conductivity (Al-Ameri and Gamadi 2019). Anderson and Fredrickson (1989) in their dynamic acid etching tests, demonstrated that higher fracture conductivity is achieved at higher temperatures. Mou et al. (2010) studied the relationship between etching pattern, conductivities and mineralogy distributions and permeability and explored the formation specification that leads to deep and narrow channels. As the result of their study, the rougher the surface the higher the conductivity is, which is because of wider fractures and deeper channels that are formed. Pournik et al. (2013) examined the effect of acid spending on fracture conductivity and fracture length. In his study, unspent acid generated a large amount of etching and partially acid system created more conductive etching pattern.
Also, fluid selection plays a crucial role on fracture conductivity. Favorable fracture conductivity can be achieved by appropriate design of acid fluid. Cash et al. (2016) performed acid conductivity experiments on samples containing a high amount of calcite. He concluded that a concentration of 15 wt% HCl has better-sustained conductivity compared with 28 wt% HCl. Also, the conductivity can be adjusted if the optimal design for acid injection time and acid concentration is performed. Pournik et al. (2007) used a different acid system for fracturing carbonate formations. Acid viscosified with surfactant, emulsified acid and acid viscosified with polymer used in their studies. There were significant differences in the fracture conductivity created between three acid systems. Their results showed that at different contact times and temperatures, conductivities from hundred to tens of thousands can be achieved. Zhang et al. (2020) performed a series of experiments on reservoir cores to study the effect of clean acid and gelled acid on carbonate formations. Rough etching created by clean acid, and channel etching created by gelled acid was observed. Accordingly, higher fracture conductivity and longer fracture length was obtained as a result of gelled acid injection.
In this study, five different fracture permeabilities ranging from 1 to 50 Darcy were used according to the published data in the experimental and field studies, and subsequently fracture conductivities between 7.5 and 1500 mD-ft were obtained.

Number of fracture stages
The horizontal well section of the sector model has a length of 3000 ft and according to this, a maximum number of 9 fracture stages (N s ) were envisaged for it. More fracture stages, due to geomechanical considerations, such as stress shadowing phenomena, seemed not to be a reasonable scenario for this system. Based on this, different cases with fracture stages of 1, 3, 5, 7 and 9 were simulated. The fractures were spaced evenly with equal spacing between them along the wellbore. To make the comparison between different cases meaningful, the location of the fractures was not changed with increasing number of stages. Only in the case of 7 stages, because of asymmetrical fracture numbers, two fracturing schemes were considered. That is, in Type-I, fractures were concentrated toward the toe and hill of the horizontal section, while in Type-II, they were concentrated toward the middle of the wellbore as shown in Fig. 5. Table 4 summarizes the full range of all fracture parameters used in this study. The simulation results demonstrated that Type-I of seven-stage fractured well performed slightly better than Type-II for all simulation cases, although the observed difference was insignificant. Therefore, only the results of the Type-I has been presented and discussed in this work.

Analysis of simulation results
A full factorial experimental design was used to generate a full combination of all fracture parameters. Based on this, a total of 600 fracturing scenarios were simulated for each reservoir model. To facilitate the simulation runs, a MATLAB computer (1) code was developed (MATLAB Software 2017b). It could automatically create the corresponding Eclipse data file to each MSFHW scenario, run the model and export the results. To evaluate the fracture performance, an improvement factor (IF) was defined, expressing the percentage gain in cumulative oil production of the fractured well compared to the non-fractured well. The corresponding formula is shown by Eq. 2.
where N pf and N pnf is the cumulative oil production of fractured and non-fractured wells, respectively.
To gain a better picture of the individual and combinational effects of fracture parameters on the well performance, practical charts were also generated. The results will be discussed in the next section.

Ranking effective parameters on well performance
Having generated a large set of simulation data (a total of 2400 simulation cases for 4 reservoir models), the response surface methodology (RSM) was employed to explore the relationship between the fracture parameters and the calculated IF. The level of impact of parameters on well performance was ranked according to the fitted model.  1 3 RSM is a useful technique to find the correlation between the interpretive variable and response variable. RSM employs different functions such as linear and full quadratic polynomial regression, to create the relationship between the output and input parameters (Al-Mudhafar and Sepehrnoori 2018;Fedutenko et al. 2012;Zubarev 2009). The full quadratic regression model with interaction terms, used in this study, has the following form: where β 0 , β i , β ij and β ii are equation coefficients and x i,j are input variables.
In this work, the input variables were fracture length, fracture conductivity and the number of fracture stages and the improvement factor (IF) was considered as the output variable. According to this, the general form of the fitted proxy model is as follows: It should be noted that the fracture width was not used in the model as a single variable, since its effect on well production has been considered together with fracture permeability by introducing the fracture conductivity into the model.

Economic analysis
Considering that MSFHW is a costly operation, optimization of fracture parameters from an economic perspective is a crucial part of this job. In other words, each of the MSFHW parameters can restrict the treatment design with respect to their operation costs. For example, fracturing cost increases considerably with increasing number of fracture stages (Rahman et al. 2014). The net present value (NPV) is a simple but very useful measure that can be used to evaluate the profitability of a project with time. In this study, considering the small size of the sector model and short production life of the reservoir, we employed a 1-year scenario for this economic exercise.
The total cost of the MSFHW job can be calculated according to Eq. 5.
where C T is the total cost of fracturing job, V a is the required volume of injected acid for two wings that is given in Eq. 6, P a is the acid price per unit volume, C Fix is a fixed cost related to transportation, crew, equipment and facilities costs and N s is the number of fracture stages. More information and details about MSFHW cost can be found in Table 5.
where x f is the fracture half-length, w f is fracture width, h f is fracture height and c L is leak-off volume coefficient.
The leak-off volume coefficient (c L ) in above formulation is considered to ensure enough volume of fracturing fluid is injected into the formation to reach ideal fracture geometry. The c L coefficient is a function of various parameters such as porosity, permeability, pressure drop across the mud cake and fluid and rock properties. Considering that leak-off volume of 1.4-3.3 times the actual fracture volume has been practiced during the common fracturing operations (Economides et al. 2002), an average value of 2 was assumed for c L coefficient in this study.
It should be mentioned that in multi-stage fracturing operations, it is usually assumed that addition of each fracture stage adds up a fixed cost of 10% to the operation cost. The cost of operation is obtained by summing up the values of items 2-4 in Table 5 that equals to 300,000$ for the case under study.
The following sequence was considered to calculate the corresponding NPV for each fracturing scenario: 1. The required volume of acid is calculated based on fracture dimensions. 2. The fixed cost of operation is calculated. (assumed 300,000 $ in this study). 3. The total cost of operation is calculated using Eq. 5. 4. The increase in cumulative oil production from fractured, compared to non-fractured, well is calculated by employing the simulation results (Eq. 7): where ΔN p,n is the increase in cumulative production after n years, N f p,n and N nf p,n are the cumulative oil production from fractured and non-fractured well after n years. 5. The revenue after acid fracturing treatment is calculated using Eq. 8. where m is the number of years evaluation is performed, and "i" is the currency escalation rate, assumed 10% in this work.

Effect of individual fracture parameters
The effect of each individual fracture parameter on well performance was studied by analyzing the trend of IF with increasing the magnitude of that parameter, whilst the other parameters were kept at their minimum, medium and maximum levels according to the ranges given in Table 4. The curves corresponding to four selected reservoir models are depicted through Figs. 6,7,8,9. Comparison of IFs between four different reservoir models shows that, first of all, the performance of fractured well is more remarkable in rocks with lower permeability. For instance, in homogenous reservoir models, when all fracture parameters are at their maximum level, IFs are about 100% and 600% for Case 1 with K = 5mD (Fig. 6) and Case 2 with k = 0.5 mD (Fig. 7), respectively. Similarly, in case of heterogeneous reservoirs, IFs are about 30% and 130% for Case 3 (Fig. 8) and Case 4 ( Fig. 9), respectively. In addition, it is noted that the multistage fractured well has considerably better performance in homogenous reservoirs in comparison with their analog heterogeneous ones, despite the average reservoir permeability is almost similar in both cases. For instance, compare Fig. 6 (Case 1) with Fig. 8 (Case 3) and Fig. 7 (Case 2) with Fig. 9 (Case 4). This, in turn, demonstrates that the distribution of rock permeability within the reservoir, especially in the regions around the fractured zones, is an important parameter affecting the fracture performance.
Analyzing the performance of fractured wells under different levels, it is observed that in all four reservoir cases the effect of any individual parameter on fracture performance is negligible when all other parameters are at their minimum level (blue curves in Fig. 6 , 7, 8, 9). As noted, the positive effect of an individual parameter becomes evident when other fracture parameters are at their medium level (orange curves). More importantly, this positive effect becomes notably more significant when other parameters are at their maximum level (green curves). In other words, there is a monotonic increase in IF as the magnitude of an individual fracture parameter increases when other parameters are at their medium and-in particular-maximum levels. For example, for the real sector model, Case 3, it Fig. 6 Effect of fracture parameters on production improvement after 1 year for case 1 1 3 is observed that when X f increases from its minimum (50 ft) to maximum (400 ft) magnitude, IF increases from 4 to 8% and from 10 to 30%, for the cases with parameters on their medium and maximum levels, respectively (Fig. 8). As another example, in Case 4, increasing the number of fracture stages from 1 to 9, results in increase in production gain from 5 to 50% when other parameters are at their medium level and from 15 to 130% when they are at their maximum Fig. 7 Effect of fracture parameters on production improvement after 1 year for case 2 Fig. 8 Effect of fracture parameters on production improvement after 1 year for case 3 level (Fig. 9). The previous observations demonstrate that the positive impact of an individual fracture parameter becomes more pronounced with increasing level of other parameters. Furthermore, there is typically a certain range for fracture parameters below which the fracturing process becomes inefficient. For instance, for the cases considered in this study, each fracture parameter showed minimal effect on well performance when all other fracture parameters were at their minimum level.
It is also interesting to note that in all four reservoir cases studied here (Figs. 6,7,8,9), increasing the fracture conductivity has the most positive effect on improving the well productivity compared to other parameters including fracture length, fracture width and number of fracture stages. This is realized by noting that the slope of corresponding line to the conductivity parameter in each reservoir case has the largest slope compared to other parameters. In other words, when fracture conductivity increases from its minimum to maximum value, IF increases to a greater extent compared to what observed for other fracture parameters.
More useful charts were also prepared from the simulation results for different reservoir cases, demonstrating the interaction effect of fracture parameters on well performance. The results are presented in Appendix A.

Level of importance of fracture parameters
To find the level of importance of each fracture parameter on well performance, the RSM technique was employed. Following this, a full quadratic proxy model was fitted to the available simulation results (IFs) for each reservoir case. Fracture length, fracture conductivity and the number of fracture stages were selected as the model variable. It should be noted that the fracture conductivity accounts for the effects of fracture width and permeability. The corresponding fitted proxy models are presented in Appendix B.
The values of the coefficient of determination of the fitted regression proxy models for each reservoir case are listed in Table 6. The estimated R 2 values verify the accuracy and reliability of the fitted models. The pareto charts corresponding to each reservoir case have also been generated and shown by Figs. 10,11,12,13. As Fig. 10 shows, in homogenous reservoir with k = 5 mD, the fracture conductivity and in the second place, the number of fracture stages are the most important parameters on enhancing the well productivity. The interaction between these two parameters also in the third place dominates the fracture performance. It is interesting to note that in this case, the impact of fracture length compared to other parameters is considerably less important.
In the second reservoir model, Case 2, which in comparison to Case 1 has a tighter rock with one order of Fig. 9 Effect of fracture parameters on production improvement after 1 year for case 4 1 3 magnitude, the number of fracture stages (N s ) dominates the well performance in the first place (Fig. 11). The fracture conductivity and then the fracture length are the next two subsequent parameters that almost with the same level of importance control the fracture performance. The observations demonstrate that in tighter reservoir systems increasing the number of fracture stages becomes significantly more important when higher well productivity is required. Furthermore, fracture length also plays more crucial role in rocks with poor flow characteristics.
Investigating the importance of fracture parameters in the heterogeneous sector model, Case 3 (Fig. 12), demonstrates more or less the same trend as observed for the homogeneous Case 1. That is, the fracture conductivity and number of fracture stages are first two important parameters, respectively, that influence the well performance. An obvious difference between two cases is the more pronounced effect of fracture length in Case 3 compared to Case 1. The importance of X f becomes even more important than that of interaction between C f and N s , observed in Case 1. This behavior can be attributed to the distribution of various rock permeabilities in Case 3, compared to a single value used in Case 1.
The results from the heterogeneous tight reservoir model (Case 4), depicted by Fig. 13, demonstrate that similar to Case 2 (low-permeability homogeneous reservoir), the number of fracture stages is the dominant parameter that controls the performance of fractured well. This is in contrast to Cases 1 and 3 (reservoirs with higher permeabilities), where the fracture conductivity was the dominating parameter. Also, the fracture length and facture conductivity are in the second and third place, respectively, affecting the well performance  . 10 The level of importance of fracture parameters in Case 1 Fig. 11 The level of importance of fracture parameters in case 2 Fig. 12 The level of importance of fracture parameters in case 3

Fig. 13
The level of importance of fracture parameters in case 4 in Case 4. The observations put more emphasis that as rock permeability decreases, increasing the number of fracture stages rather than fracture conductivity become more important and profitable to expedite the well production rate.

Optimization under economic considerations
In the previous sections, the performance of fractured well was investigated from a technical point of view. In other words, maximizing the well productivity by changing the fracture parameters was considered as the main object of fracturing job. As discussed earlier, due to excessive cost of multi-stage fracturing operations, it is vital to optimize the fracture design from an economic perspective. Accordingly, to find the optimum economic fracture parameters, maximizing NPV was used as the final criteria to select the best scenario. The evaluations were performed for all four reservoir models based on the procedures discussed previously in Sect. 3.4. As mentioned earlier considering the short production life cycle of the reservoir, a 1-year production scheme was used for the NPV calculations. Table 7 shows the results of the economic exercise performed, where the optimum fracture parameters and the corresponding maximum NPV for each reservoir model have been determined.
It should be mentioned that when only the technical considerations were important, i.e. evaluations presented in Sect. Effect of Individual Fracture Parameters, the maximum MSFHW productivity for all four reservoir models was attained when all fracture parameters were at their maximum levels. That is, X f = 400 ft, W f = 0.36 in, C f = 1500 md-ft and N s = 9. However, in case of economic considerations, the optimum fracture parameters for some cases lie in an obviously different range. For instance, the optimum fracture length for the Case 3 is 50 ft and optimum fracture width and conductivity for Case 4 are 0.09 in and 375 md-ft, respectively (Table 7). Such observations place emphasis on the significance of careful economic analysis within reasonable evaluation time intervals, where an optimum fracture design is demanded.
To provide a better picture of importance of economic analysis for fracture designs, the well production IFs corresponding to technical (maximum) and economic (optimum) fracture designs have been compared in Table 8. It is interesting to note that, especially for the heterogeneous reservoir cases (Cases 3 and 4), there is an obvious difference between IFs of optimum and maximum design scenarios. It is observed that for the real reservoir model (base Case 3), the maximum IFs can reach about 29%, but the economic design restricts this value to about 11%. The observations underline the strict economic analyses required to achieve the most profitable fracture design for field applications.

Conclusions
The work presented in this paper focused on optimization of fracture parameters in MSFHWs in low-permeability heavyoil systems. Two types of homogeneous and heterogeneous sector reservoir models with permeabilities ranging between 0.5 and 5 mD were used and about 2400 simulation scenarios were run. The RSM methodology was then employed to rank the level of importance of fracture parameters. The optimum fracture parameters were also determined by conducting an economic exercise by calculating NPVs corresponding to each fracture design. The main conclusions from this study are as follows: • For the reservoir cases with average permeability of 5 mD, fracture conductivity and number of fracture stages, were in turn the main effective parameters improving the well productivity. As rock permeability decreased, i.e. 0.5 mD, the number of fracture stages became the dominant parameter controlling the fractured well performance. • In tight reservoir systems, increasing the fracture length became a beneficial practice to reach higher production rates. However, it should be mentioned that if the fracture conductivity is low, it is a good practice to keep the fracture length shorter to avoid additional costs. • Considering reservoir heterogeneity, especially distribution of rock permeability in the stimulated area around the wellbore, showed a significant role in determining the optimum fracture design. For the tight heterogeneous sector model, the fracture length had a more posi- tive effect on well productivity compared to fracture conductivity, whereas the later parameter was always the dominating factor over the former one in homogeneous systems. • The economic analysis of the fracturing operation played a key role in designing the optimum fracture parameters. More specifically determination of maximum fracture length and number of fracture stages are crucially affected by profitability of the operation.
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/.

Appendix A
The "main effects plot" examine differences between level means (average input values) for various parameters. These graphs are prepared based on analysis of variance (ANOVA). The response means (average output values) for each level of parameter are plotted and connected by a line. In these charts, the value of each parameter increases while the value of the other parameters are kept at average values. The baseline represents the value of the response level in the mean value of all parameters. Also, a reference line (dashed line) is drawn at the overall mean. The reference line represents the value of the response level in the mean value of all parameters. These charts verified our results through Effect of Individual Fracture Parameters and Level of importance of fracture parameters sections (Figs. 14,15,16,17).
The following contour plots show the interaction effect of fracture parameters on well performance (IF) for different reservoir cases. As the highlighted region becomes darker, the effect of parameters become more favorable and IF increases (Figs. 18,19,20,21).