Design and Performance Analysis of Dry Gas Fishbone Wells for Lower Carbon Footprint

: Multilateral well drilling technology has recently assisted the drilling industry in improving borehole contact area and reducing operation time, while maintaining a competitive cost. The most advanced multilateral well drilling method is Fishbone drilling (FbD). This method has been utilized in several hydrocarbon ﬁelds worldwide, resulting in high recovery enhancement and reduced carbon emissions from drilling. FbD involves drilling several branches from laterals and can be considered as an alternative method to hydraulic fracturing to increase the stimulated reservoir volume. However, the expected productivity of applying a Fishbone well from one ﬁeld to another can vary due to various challenges such as Fishbone well design, reservoir lithology, and accessibility. Another challenge is the lack of existing analytical models and the effect of each Fishbone parameter on the cumulative production, as well as the interaction between them. In this paper, analytical and empirical productivity models were modiﬁed for FbD in a dry gas reservoir. The modiﬁed analytical model showed a higher accuracy with respect to the existing model. It was also compared with the modiﬁed empirical model, which proved its higher accuracy. Finally, machine learning algorithms were developed to predict FbD productivity, which showed close results with both analytical and empirical models.


Introduction
A multilateral well configuration encompasses the deployment of multiple horizontal laterals branching off from a singular vertical wellbore, with the objective of enhancing production efficiency [1]. The utilization of multilateral wells aims to increase the drainage area and thus augment the production rate [2]. Although the drilling cost associated with multilateral wells may be higher, this cost is offset by the improved reservoir recovery resulting from enhanced wellbore-reservoir contact. There has been extensive research directed towards optimizing multilateral drilling costs and designs [3,4]. One of the latest advancements in multilateral drilling technology is Fishbone Drilling (FbD) [5]. FbD is comprised of a series of micro-branches drilled in divergent directions from a primary horizontal wellbore, with a designated inter-branch spacing as shown in Figure 1 [6]. The implementation of FbD can maximize reservoir contact area and facilitate access to previously unreachable geological formations, thereby augmenting cumulative production [7]. FbD has been shown to have higher productivity compared to hydraulic fracturing in relatively low permeability formations [8]. FbD is a desirable solution in the light of economic, environmental, and regulatory limitations [9]. A recent study conducted by THREE60 Energy found that Fishbone drilling methods result in significantly lower CO2 emissions The optimization of Fishbone well parameters has been a subject of investigation various researchers, who have sought to comprehend the effect on productivity and mulative production [14]. This includes the impact of branch number, branch length ter-branch distance, and the angle between the branch and main lateral borehole. Fishbone configuration has been implemented in the Daqing Peripheral oilfield [8] an Middle East oilfield [14]. Numerical simulations demonstrate that well productivity is hanced with an increase in branch number, reaching an optimum point. The optim branch number represents the least number of branches required to achieve maxim productivity increase. Analogous findings were observed with increasing branch len with an optimum length of 300 m, beyond which productivity remained constant. optimal branch angle with respect to the main lateral was found to be 20° and 30° in different cases, influenced by the control area and its relation to angle changes. The in branch distance is still under investigation, as it pertains to multiple factors such as re voir thickness, petrophysical properties, drilling Kickoff Point (KOP) selection, and bu up severity challenges. To maximize reservoir volume, the drainage area should b sufficient width to eliminate interference between adjacent branches.
Successful applications of Fishbone well configurations have been reported in sev case studies across the globe. In Venezuela, the Zuata field, which comprises non-ho geneous sand, has been documented as the first successful application of Fishbone des (FbD) [15]. The first FbD application in Russia was reported in the Nizhnekhet 1 fiel very thin formation [16]. Four additional case studies from Russia have been documen in the literature, including sandstone formations in the Vostochno-Messoyakskoye shore field [17], low permeable formations in the East Messoyakha field [18], the Van oilfield in Krasnoyarsk, Eastern Siberia [19], and the Srednebotuobinskoye oil and condensate field in Eastern Siberia [20]. In the UAE, successful FbD case studies have b reported in naturally fractured reservoirs [21,22], and an innovative Fishbone infill w was applied in the Canadian oil sands play [23]. Research into applications of Fishb well configurations in areas such as geothermal drilling, enhanced oil recovery, ste assisted gravity drainage, and underground saltwater disposal are ongoing [24,25]. The optimization of Fishbone well parameters has been a subject of investigation by various researchers, who have sought to comprehend the effect on productivity and cumulative production [14]. This includes the impact of branch number, branch length, inter-branch distance, and the angle between the branch and main lateral borehole. The Fishbone configuration has been implemented in the Daqing Peripheral oilfield [8] and a Middle East oilfield [14]. Numerical simulations demonstrate that well productivity is enhanced with an increase in branch number, reaching an optimum point. The optimum branch number represents the least number of branches required to achieve maximum productivity increase. Analogous findings were observed with increasing branch length, with an optimum length of 300 m, beyond which productivity remained constant. The optimal branch angle with respect to the main lateral was found to be 20 • and 30 • in two different cases, influenced by the control area and its relation to angle changes. The inter-branch distance is still under investigation, as it pertains to multiple factors such as reservoir thickness, petrophysical properties, drilling Kickoff Point (KOP) selection, and build-up severity challenges. To maximize reservoir volume, the drainage area should be of sufficient width to eliminate interference between adjacent branches.
Successful applications of Fishbone well configurations have been reported in several case studies across the globe. In Venezuela, the Zuata field, which comprises nonhomogeneous sand, has been documented as the first successful application of Fishbone design (FbD) [15]. The first FbD application in Russia was reported in the Nizhnekhet 1 field, a very thin formation [16]. Four additional case studies from Russia have been documented in the literature, including sandstone formations in the Vostochno-Messoyakskoye onshore field [17], low permeable formations in the East Messoyakha field [18], the Vankor oilfield in Krasnoyarsk, Eastern Siberia [19], and the Srednebotuobinskoye oil and gas condensate field in Eastern Siberia [20]. In the UAE, successful FbD case studies have been reported in naturally fractured reservoirs [21,22], and an innovative Fishbone infill well was applied in the Canadian oil sands play [23]. Research into applications of Fishbone well configurations in areas such as geothermal drilling, enhanced oil recovery, steam-assisted gravity drainage, and underground saltwater disposal are ongoing [24,25].
The production rate of oil and gas from a reservoir can be predicted through reservoir deliverability analysis. The bottom-hole pressure has a significant effect on the production rate, regardless of the completion type and artificial lift methods employed. Analytical models for well deliverability exist to study the relationship between the production rate and bottom-hole pressure [26]. These models are dependent on various factors such as reservoir pressure, boundary type, distance, permeability (vertical, horizontal, and relative permeability), wellbore radius, near-wellbore effect, thickness of the pay zone, and reservoir fluid properties [27]. The prediction of productivity is also dependent on well geometry (e.g., vertical, horizontal, multilateral) and the interrelationship between the parameters mentioned is known as the Inflow Performance Relationship (IPR).
This paper presents advanced approaches to quantify the production rates of some Fishbone well configurations. Analytical, empirical, and data-driven approaches were adopted to develop productivity prediction models in a dry gas reservoir with good petrophysical properties. The developed models consider the effect of the number and length of the branches, the distance between two adjacent branches, the inclination angle, and the permeability anisotropy. The modified analytical model was based on modifying the existing models and some geometrical assumptions. The empirical model was modified using the Grey Wolf Algorithm optimization, and the data-driven models were based on the artificial neural network and support vector regression methods optimized with an advanced genetic algorithm.

Background
The implementation of Fishbone well design (FbD) presents numerous difficulties in terms of wellbore analysis and production flow studies due to its complex combination of boreholes with differing well geometries [28]. These difficulties result in the development of analytical models based on various assumptions, leading to potential inaccuracies in production rate predictions. In an effort to mitigate such errors, numerical simulations are often employed, however, these too are subject to complications. As the progress of FbD continues to advance, the development of new models and approaches for the analysis of productivity in Fishbone wells becomes increasingly crucial.
In regard to the analysis of production flow rate, two different approaches exist, analytical and empirical. Analytical approaches are derived from the solution of Darcy's equation for different flow regimes, while empirical approaches are based on experimental work to develop correlations between drawdown and flow rate [29].
The Inflow Performance Relationship (IPR) for Fishbone wells is challenging to model due to the complexity of its geometries. Two different flow regions exist, with one representing the drilled region where the flow between branches is assumed to be pseudosteady-state and pseudo-linear flow, and the other representing the formation matrix where the flow is assumed to be pseudo-steady-state radial flow. The inner region denotes the drained volume, while the outer region symbolizes the non-stimulated volume [30] (See Appendix A).
The deliverability of a uniformly distributed number of holes in the inner region for a gas reservoir was presented by [30] (pp. . The model is stated as: Equation (1) is used for the inner regions where the formation is drilled using Fishbone branches; Equation (2) is used for the outer region: In the above equations: L i length of the branch (ft), n number of the branches, r w i radius of the branch (ft), y b i distance between adjacent branches (ft) (These parameters represent the Fishbone well design configuration). s i Skin factor, µ g : Gas viscosity (cp), z Gas compressibility factor, T Reservoir temperature ( • R), p i Initial reservoir pressure (psia), p w f Flowing bottom-hole pressure r = rw (psia), r w Well radius (ft), s Skin factor, r e Reservoir boundary (ft), p r Wellbore pressure (psi), h Pay thickness (ft), q g Gas flow rate at pressure p (bbl/day), D Non-darcy flow coefficient (d/Mscf), C A Stabilized performance coefficient (MMScf/D.psia 2m or MMScf/D.(psia 2 /cp) m ), m Dimensionless deliverability exponent, defined as the line's inverse slope on a log-log plot of the change in pressuresquared or pseudo-pressure versus gas flow rate, p pl Average pressure at the edge of the inner region (psi), γ Fluid specific gravity, I ani is the relationship between the horizontal and vertical permeability, which is presented as: where: k H Horizontal permeability (mD), k V Vertical permeability (mD), r pl Equivalent radius of the inner region, which can be calculated using the following equation: Due to the complexity of solving analytical models, caused mainly by the lack of accurate PVT data, IPR correlations were introduced to the petroleum industry. Ref. [31] generated an empirical correlation to model IPR in a Fishbone type multilateral dry gas well. Their model was validated using a commercial well performance software simulation with real field data for production history matching. The model does not take into consideration the effect of the distance between the adjacent branches on productivity. This model does not yield accurate results when the number of branches increases.
Several papers and research studies report on the productivity prediction of vertical and horizontal wells using different machine learning algorithms [32,33]. It has been recently used as a powerful tool to predict the well performance in complex formations [34]. Ref. [35] developed algorithms based on the Neural-based Decision Tree (NDT) learning model. A higher-order neural network was applied by [36], and a surrogate model was suggested by [37]. The other models are based on integrating supervised machine learning with a well-calibrated bias [38], and experimental data and/or simulation results [39]. Ref. [40] developed a machine learning approach for shale gas horizontal well productivity.
Ref. [41] proposed three methods for artificial intelligence as a new approach to investigate the productivity in Fishbone well, the inputs are the flowing bottom hole pressure (P wf ), the formation permeability ratio (k H /k V ), and the length of each lateral (Li). The model output was the flow rate fraction. This study is limited in terms of the other Fishbone parameters such as the number of branches and distance between them. The three applied algorithms were Artificial Neural Network (ANN) model, with 97% accuracy, the Adaptive Neuro-Fuzzy Inference System (ANFIS) with 98% R score, and the Radial Basis Function (RBF) network with a 98% coefficient of correlation.

Simulation Model Description
Numerical simulation plays a crucial role in the analysis of well deliverability when field data is scarce or limited, which is the case for a newly developed technology such as fishbone drilling. For this reason, in our study, we consider the results from the numerical simulation to be the benchmark for comparison and validation of our models.
A commercial well performance software was used to evaluate the performance of Fishbone wells at different operating conditions. More than 250 different cases were examined. Actual wellbore conditions and real reservoir data were utilized to perform the study. The reservoir model was examined and validated by an experienced production engineer. The reservoir model consists of 62 × 21 × 11 grid cells, while a 3D Cartesian grid system was used. The Fishbone well was drilled at the center of a gas reservoir that has dimensions of 20,000 × 10,000 × 750 feet. The Fishbone well performance was studied for around 20 years of production, representing a typical field life span in the area. Further descriptions about the reservoir models have been reported previously by [31,41]. In all simulation runs, the well production and bottomhole pressure profiles were recorded. In addition, wide ranges of permeability ratio (k H /k V ), the distance between adjacent branches, and the number of branches were used. Table 1 lists the statistical parameters for all reservoir models studied in this work. The production rate varied between 197,900 and 0 SCF/DAY, representing the Absolute Open Flow (AOF) and well shut-in conditions, respectively. A permeability anisotropy factor (k H /k V ) between 1000 and 1 was used, which represents very tight vertical permeability and equivalent horizontal and vertical permeability (uniform and homogenous reservoir), respectively. Moreover, flowing well pressure between 14.7 psi and 4800 psi was used, indicating AOF and well shut-in conditions.

Model Elaboration
The results generated from the numerical simulations are presented in this section and are used as the input for the analytical, empirical, and data-driven model elaboration. The numerical simulation results are compared with the existing analytical model which is modified and compared with the modified empirical model.

Analytical IPR Model for Fishbone Well
Based on the literature review presented previously, we modified the analytical models for FbD productivity prediction following the work of [30] (pp. 37-81). The aim was to apply the existing model to the reservoir simulation results and investigate the applicability and reliability of the numerical simulation outputs. The comparison is based on each parameter change, and the output is presented in a chart format.

. Number of Branches
To evaluate the impact of the number of branches on cumulative production, a constant permeability anisotropy of 10 was established while maintaining a fixed length of 3100 feet and an interbranch distance of 2600 feet. The only varying parameter was the number of branches. The results depicted in Figure 2 indicate an exponential increase in flow rates obtained from the analytical model in comparison to a slight increase observed in the numerical simulation results with the increase in the number of branches. These findings verify the limitations of the existing analytical model for a higher number of branches, specifically above four branches.
Fuels 2023, 4, FOR PEER REVIEW 6 applicability and reliability of the numerical simulation outputs. The comparison is based on each parameter change, and the output is presented in a chart format.

Number of Branches
To evaluate the impact of the number of branches on cumulative production, a constant permeability anisotropy of 10 was established while maintaining a fixed length of 3100 feet and an interbranch distance of 2600 feet. The only varying parameter was the number of branches. The results depicted in Figure 2 indicate an exponential increase in flow rates obtained from the analytical model in comparison to a slight increase observed in the numerical simulation results with the increase in the number of branches. These findings verify the limitations of the existing analytical model for a higher number of branches, specifically above four branches.

Length of Branches
In our analysis of the influence of branch length on cumulative production, a permeability anisotropy of 10 was maintained while keeping the number of branches constant at 6 and inter-branch distance fixed at 2600 feet. The sole variable in this scenario was the length of the branches. Our numerical simulation and the existing analytical model both depict an increase in well productivity with an increase in branch length. However, the gradients of the productivity growth chart for both models remain unchanged. The existing analytical model demonstrates an overestimation of the impact of the length of branches on flow rate as illustrated in Figure 3, when compared to the results obtained from the numerical simulation.

Permeability Anisotropy
In this study, we evaluated several scenarios where the length of the branches was fixed at 3100 feet, the inter-branch distance was maintained at 2600 feet, and the number

Length of Branches
In our analysis of the influence of branch length on cumulative production, a permeability anisotropy of 10 was maintained while keeping the number of branches constant at 6 and inter-branch distance fixed at 2600 feet. The sole variable in this scenario was the length of the branches. Our numerical simulation and the existing analytical model both depict an increase in well productivity with an increase in branch length. However, the gradients of the productivity growth chart for both models remain unchanged. The existing analytical model demonstrates an overestimation of the impact of the length of branches on flow rate as illustrated in Figure 3, when compared to the results obtained from the numerical simulation. applicability and reliability of the numerical simulation outputs. The comparison is based on each parameter change, and the output is presented in a chart format.

Number of Branches
To evaluate the impact of the number of branches on cumulative production, a constant permeability anisotropy of 10 was established while maintaining a fixed length of 3100 feet and an interbranch distance of 2600 feet. The only varying parameter was the number of branches. The results depicted in Figure 2 indicate an exponential increase in flow rates obtained from the analytical model in comparison to a slight increase observed in the numerical simulation results with the increase in the number of branches. These findings verify the limitations of the existing analytical model for a higher number of branches, specifically above four branches.

Length of Branches
In our analysis of the influence of branch length on cumulative production, a permeability anisotropy of 10 was maintained while keeping the number of branches constant at 6 and inter-branch distance fixed at 2600 feet. The sole variable in this scenario was the length of the branches. Our numerical simulation and the existing analytical model both depict an increase in well productivity with an increase in branch length. However, the gradients of the productivity growth chart for both models remain unchanged. The existing analytical model demonstrates an overestimation of the impact of the length of branches on flow rate as illustrated in Figure 3, when compared to the results obtained from the numerical simulation.

Permeability Anisotropy
In this study, we evaluated several scenarios where the length of the branches was fixed at 3100 feet, the inter-branch distance was maintained at 2600 feet, and the number

Permeability Anisotropy
In this study, we evaluated several scenarios where the length of the branches was fixed at 3100 feet, the inter-branch distance was maintained at 2600 feet, and the number of branches was held constant at 6, while varying the permeability anisotropy. The results displayed in Figure 4 indicate a correlation between increased reservoir anisotropy and enhanced well productivity. The outputs from both the existing analytical models and Fuels 2023, 4 98 the numerical simulations are highly congruent, with a common value at the intersection between the two models represented by a permeability anisotropy I ani of 40-50.
Fuels 2023, 4, FOR PEER REVIEW 7 of branches was held constant at 6, while varying the permeability anisotropy. The results displayed in Figure 4 indicate a correlation between increased reservoir anisotropy and enhanced well productivity. The outputs from both the existing analytical models and the numerical simulations are highly congruent, with a common value at the intersection between the two models represented by a permeability anisotropy Iani of 40-50.

Distance between Adjacent Branches
In this examination, we analyzed the impact of inter-branch distance on cumulative production by maintaining a constant number of branches (6), permeability anisotropy (10), and length of branches (3100 feet). The results obtained from the existing analytical model and numerical simulations displayed stark disparities as demonstrated in Figure  5. The trend observed in the numerical simulation depicts a slight decrease, whereas the analytical model exhibits an increase. These findings affirm that the existing analytical model is not suitable for studying the effect of inter-branch distance variations. The development of the analytical model assumes that as inter-branch distance increases, the drainage area expands. This theory, however, is limited when the inter-branch distance is less than double the drainage radius. An increase in inter-branch distance results in an increase in cumulative production due to the enlarged drainage area. When the inter-branch distance equals double the drainage radius, the cumulative production remains constant and experiences a slight decrease in regions with poor reservoir quality properties. The simulation is configured with inter-branch distances greater than double the drainage radius to prevent any interference with the drainage volume of each branch. As such, the modification of the existing model is based on two parameters:

Distance between Adjacent Branches
In this examination, we analyzed the impact of inter-branch distance on cumulative production by maintaining a constant number of branches (6), permeability anisotropy (10), and length of branches (3100 feet). The results obtained from the existing analytical model and numerical simulations displayed stark disparities as demonstrated in Figure 5. The trend observed in the numerical simulation depicts a slight decrease, whereas the analytical model exhibits an increase. These findings affirm that the existing analytical model is not suitable for studying the effect of inter-branch distance variations.
Fuels 2023, 4, FOR PEER REVIEW 7 of branches was held constant at 6, while varying the permeability anisotropy. The results displayed in Figure 4 indicate a correlation between increased reservoir anisotropy and enhanced well productivity. The outputs from both the existing analytical models and the numerical simulations are highly congruent, with a common value at the intersection between the two models represented by a permeability anisotropy Iani of 40-50.

Distance between Adjacent Branches
In this examination, we analyzed the impact of inter-branch distance on cumulative production by maintaining a constant number of branches (6), permeability anisotropy (10), and length of branches (3100 feet). The results obtained from the existing analytical model and numerical simulations displayed stark disparities as demonstrated in Figure  5. The trend observed in the numerical simulation depicts a slight decrease, whereas the analytical model exhibits an increase. These findings affirm that the existing analytical model is not suitable for studying the effect of inter-branch distance variations. The development of the analytical model assumes that as inter-branch distance increases, the drainage area expands. This theory, however, is limited when the inter-branch distance is less than double the drainage radius. An increase in inter-branch distance results in an increase in cumulative production due to the enlarged drainage area. When the inter-branch distance equals double the drainage radius, the cumulative production remains constant and experiences a slight decrease in regions with poor reservoir quality properties. The simulation is configured with inter-branch distances greater than double the drainage radius to prevent any interference with the drainage volume of each branch. As such, the modification of the existing model is based on two parameters: • Drainage surface: The development of the analytical model assumes that as inter-branch distance increases, the drainage area expands. This theory, however, is limited when the inter-branch distance is less than double the drainage radius. An increase in inter-branch distance results in an increase in cumulative production due to the enlarged drainage area. When the inter-branch distance equals double the drainage radius, the cumulative production remains constant and experiences a slight decrease in regions with poor reservoir quality properties. The simulation is configured with inter-branch distances greater than double the drainage radius to prevent any interference with the drainage volume of each branch. As such, the modification of the existing model is based on two parameters:

•
Drainage surface: The drainage surface is a part of the stimulated radius "r pl ". The r pl of the existing correlation considers that by increasing the distance between adjacent branches, the drainage surface increases. To develop that, we take an example of a Fishbone well with five branches, as shown in Figure 6.
The drainage surface is a part of the stimulated radius "rpl". The rpl of the existing correlation considers that by increasing the distance between adjacent branches, the drainage surface increases. To develop that, we take an example of a Fishbone well with five branches, as shown in Figure 6.
= √ ( + 1)2 The developed stimulated radius assumes that the drainage volume is related to the drainage radius and not the distance between adjacent branches. In the case where the double drainage radius is greater than the distance between adjacent branches, Equation (13) can be used instead of Equation (7).

•
High-pressure reservoir effect: We noticed that the existing analytical models do not consider the effect of high reservoir pressure. For gas reservoirs, when the pressure exceeds 3000 psi, the pressure function (1/μg Bg) becomes nearly constant, and thus can be taken outside the integral. This function will replace (1/μg z) used in low pressure reservoirs [42]. This approximation is called pressure-approximation method. Equation (1) is modified accordingly First, the square function needs to be deleted from the pressure difference (this represents the gas pseudo-pressure integration): Second, the constant pressure function (1/μg Bg) is introduced into the equation (adjustments required): The development of r pl is as follows: The developed stimulated radius assumes that the drainage volume is related to the drainage radius and not the distance between adjacent branches. In the case where the double drainage radius is greater than the distance between adjacent branches, Equation (13) can be used instead of Equation (7).

•
High-pressure reservoir effect: We noticed that the existing analytical models do not consider the effect of high reservoir pressure. For gas reservoirs, when the pressure exceeds 3000 psi, the pressure function (1/µ g B g ) becomes nearly constant, and thus can be taken outside the integral. This function will replace (1/µ g z) used in low pressure reservoirs [42]. This approximation is called pressure-approximation method. Equation (1) is modified accordingly First, the square function needs to be deleted from the pressure difference (this represents the gas pseudo-pressure integration): Second, the constant pressure function (1/µ g B g ) is introduced into the equation (adjustments required): In order to apply the p-approximation of gas deliverability equations, the gas volume factor B g should be introduced, and the gas viscosity needs to be calculated. The correlations are shown in Appendix B.

IPR Empirical Correlation for Fishbone Well
This section aims to generate a new IPR empirical correlation for a Fishbone well for a conventional dry gas reservoir. The non-linearity between the parameters needs to be solved by either genetic algorithms or evolutionary algorithms. The multifeatured non-linear regression is solved using the Grey Wolf Algorithm (GWA) algorithm which was first introduced in 2014 [43]. The algorithm has been widely used in the petroleum industry [44]. The grey wolf's behavior mainly inspired the algorithm. In nature, these wolves live in packs of 5-12 wolves which have a hierarchy of dominance. The wolf that governs the pack is called alpha (α), and they are responsible for deciding on the hunt. The whole pack must follow these decisions. Beta (β) wolves come next in the hierarchy, and these are subordinate wolves helping the alpha make decisions. The beta wolves reinforce the alpha decisions through the pack. The next step is to find the delta (δ) wolves. They are subordinate of alphas and betas, but still, dominate omegas ( In order to apply the p-approximation of gas deliverability equations, the gas volume factor Bg should be introduced, and the gas viscosity needs to be calculated. The correlations are shown in Appendix B.

IPR Empirical Correlation for Fishbone Well
This section aims to generate a new IPR empirical correlation for a Fishbone well for a conventional dry gas reservoir. The non-linearity between the parameters needs to be solved by either genetic algorithms or evolutionary algorithms. The multifeatured nonlinear regression is solved using the Grey Wolf Algorithm (GWA) algorithm which was first introduced in 2014 [43]. The algorithm has been widely used in the petroleum industry [44]. The grey wolf's behavior mainly inspired the algorithm. In nature, these wolves live in packs of 5-12 wolves which have a hierarchy of dominance. The wolf that governs the pack is called alpha (α), and they are responsible for deciding on the hunt. The whole pack must follow these decisions. Beta (β) wolves come next in the hierarchy, and these are subordinate wolves helping the alpha make decisions. The beta wolves reinforce the alpha decisions through the pack. The next step is to find the delta (δ) wolves. They are subordinate of alphas and betas, but still, dominate omegas (ϒ) which have the lowest level in the hierarchy [43]. This hierarchy and the hunting strategies mentioned in [43] have been modeled mathematically to produce the Grey Wolf Algorithm. First, a random set of solutions is generated to represent the location of these wolves. Each of these solutions is evaluated using the objective function-root mean square error for our case. Second, the solutions are ranked, and the three fittest solutions are assigned as positions to the wolves' alpha (α), beta (β), and delta (δ), while the rest of the solutions are assigned to the omegas (ϒ). The wolves are looking for prey in the hunting process, which is the optimum solution in the mathematical model.
The modified correlation will be as follows: a, b, c, d, e, g, and h are the parameters that calibrate our model to ensure the accuracy between the inputs (dimensionless BHP, anisotropy, distance between adjacent branches, and the length of the branch) and the output (dimensionless flow rate). These parameters are listed in Table 2. The GWA algorithm was used to find the optimum values of the coefficient. The coefficient of correlation (R-Score) is 98%, as shown in Figure 7: ) which have the lowest level in the hierarchy [43]. This hierarchy and the hunting strategies mentioned in [43] have been modeled mathematically to produce the Grey Wolf Algorithm. First, a random set of solutions is generated to represent the location of these wolves. Each of these solutions is evaluated using the objective function-root mean square error for our case. Second, the solutions are ranked, and the three fittest solutions are assigned as positions to the wolves' alpha (α), beta (β), and delta (δ), while the rest of the solutions are assigned to the omegas ( p-approximation of gas deliverability equations, the gas volume uced, and the gas viscosity needs to be calculated. The correladix B. ion for Fishbone Well generate a new IPR empirical correlation for a Fishbone well for eservoir. The non-linearity between the parameters needs to be algorithms or evolutionary algorithms. The multifeatured nond using the Grey Wolf Algorithm (GWA) algorithm which was 3]. The algorithm has been widely used in the petroleum indusehavior mainly inspired the algorithm. In nature, these wolves es which have a hierarchy of dominance. The wolf that governs ), and they are responsible for deciding on the hunt. The whole ecisions. Beta (β) wolves come next in the hierarchy, and these elping the alpha make decisions. The beta wolves reinforce the the pack. The next step is to find the delta (δ) wolves. They are d betas, but still, dominate omegas (ϒ) which have the lowest ]. This hierarchy and the hunting strategies mentioned in [43] ematically to produce the Grey Wolf Algorithm. First, a random ed to represent the location of these wolves. Each of these solue objective function-root mean square error for our case. Second, , and the three fittest solutions are assigned as positions to the ), and delta (δ), while the rest of the solutions are assigned to the re looking for prey in the hunting process, which is the optimum ical model. tion will be as follows: e the parameters that calibrate our model to ensure the accuracy nsionless BHP, anisotropy, distance between adjacent branches, ch) and the output (dimensionless flow rate). These parameters the modified empirical model. was used to find the optimum values of the coefficient. The co--Score) is 98%, as shown in Figure 7: ). The wolves are looking for prey in the hunting process, which is the optimum solution in the mathematical model.
The modified correlation will be as follows: a, b, c, d, e, g, and h are the parameters that calibrate our model to ensure the accuracy between the inputs (dimensionless BHP, anisotropy, distance between adjacent branches, and the length of the branch) and the output (dimensionless flow rate). These parameters are listed in Table 2. The GWA algorithm was used to find the optimum values of the coefficient. The coefficient of correlation (R-Score) is 98%, as shown in Figure 7: The coefficient of correlation represented in Figure 7 demonstrates that the variable is highly correlated, so the modified correlation is statistically valid. To ensure that the correlation is physically valid, the signs (negative or positive) of the coefficient are related to the effect of the Fishbone well configuration parameters on the flow rate increase or decrease, which was added to the developed algorithm with possible intervals for each coefficient to account for the physical significance of the problem.
To validate the results, a comparison between the numerical simulation, modified empirical correlation, and modified analytical model flow rates was conducted. This was carried out by varying the number of the branches, length of the branches, anisotropy, and the distance between these branches. The coefficient of correlation represented in Figure 7 demonstrates that the variable is highly correlated, so the modified correlation is statistically valid. To ensure that the correlation is physically valid, the signs (negative or positive) of the coefficient are related to the effect of the Fishbone well configuration parameters on the flow rate increase or decrease, which was added to the developed algorithm with possible intervals for each coefficient to account for the physical significance of the problem.
To validate the results, a comparison between the numerical simulation, modified empirical correlation, and modified analytical model flow rates was conducted. This was carried out by varying the number of the branches, length of the branches, anisotropy, and the distance between these branches.

IPR Data-Driven Models for Fishbone Well
The robust data-driven models are used to predict the flow rate based on the numerical simulation results and the inputs include anisotropy (Iani), the number of the branches (n), the distance between adjacent branches (ybi), the length of each branch (Li), and the BHP pressure divided by the reservoir pressure.
The dataset is obtained from the same numerical simulation used in this study; two powerful algorithms were developed for the productivity prediction, the Support Vector Regression optimized by the Genetic Algorithm (SVR-GA) and the Artificial Neuron Network (ANN).

Artificial Neural Network (ANN)
ANN modeling can mimic human brain behavior and predict the output of a complex function from the inputs through learning [45,46]. Ref. [47] introduced the composition of the ANNs as a learning algorithm presented as an architecture of layers composed of neurons that have transfer functions embedded in them and connected by weights until the signal arrives at the predicted output [48].
The resulting regression for the ANN showed a high accuracy of 99.05% as shown in Figure 8, with a root mean square error of 0.0205, and the mean absolute error of 0.033 (see Figure 9).

IPR Data-Driven Models for Fishbone Well
The robust data-driven models are used to predict the flow rate based on the numerical simulation results and the inputs include anisotropy (I ani ), the number of the branches (n), the distance between adjacent branches (y bi ), the length of each branch (L i ), and the BHP pressure divided by the reservoir pressure.
The dataset is obtained from the same numerical simulation used in this study; two powerful algorithms were developed for the productivity prediction, the Support Vector Regression optimized by the Genetic Algorithm (SVR-GA) and the Artificial Neuron Network (ANN).

Artificial Neural Network (ANN)
ANN modeling can mimic human brain behavior and predict the output of a complex function from the inputs through learning [45,46]. Ref. [47] introduced the composition of the ANNs as a learning algorithm presented as an architecture of layers composed of neurons that have transfer functions embedded in them and connected by weights until the signal arrives at the predicted output [48].
The resulting regression for the ANN showed a high accuracy of 99.05% as shown in Figure 8, with a root mean square error of 0.0205, and the mean absolute error of 0.033 (see Figure 9).

SVR Optimized by the Genetic Algorithm (SVR-GA)
Ref. [49] stated that Support Vector Regression (SVR) is a subset of the Support Vector Machine, which is a machine learning algorithm that performs regression with high efficiency. It was developed by [50]. SVR is different from ordinary regression, where it minimizes the generalized error bound that gathers the training error and a regularization term instead of minimizing the error between the output and the original values [51].
Since the first introduction of Genetic Algorithms (GA) in the sixties to solve both linear and convex optimization along with non-linear simultaneous equations by [52], Evolutionary Algorithms (EA) have played an essential role in advancing numerical solutions and computations. The most famous and widely used algorithm among EA is the Genetic Algorithm which was first introduced by [53]. The main mechanism behind GA is the ability to mimic biological evolution. GA is a population-based algorithm, meaning that it creates an initial set of random solutions called "Chromosomes" composed of genes as the elements of the solution; then, these solutions are subject to a "Fitness Function" defined by the user and, in our case is the performance of the SVR. To obtain the best solution for the optimization, GA used what are called genetic operators [54], that keep interacting with the solution until it converges toward the global optimum.
The SVR-GR model exhibits superior accuracy compared to the ANN model due to the utilization of a robust Genetic Algorithm, effectively reducing the error between the

SVR Optimized by the Genetic Algorithm (SVR-GA)
Ref. [49] stated that Support Vector Regression (SVR) is a subset of the Support Vector Machine, which is a machine learning algorithm that performs regression with high efficiency. It was developed by [50]. SVR is different from ordinary regression, where it minimizes the generalized error bound that gathers the training error and a regularization term instead of minimizing the error between the output and the original values [51].
Since the first introduction of Genetic Algorithms (GA) in the sixties to solve both linear and convex optimization along with non-linear simultaneous equations by [52], Evolutionary Algorithms (EA) have played an essential role in advancing numerical solutions and computations. The most famous and widely used algorithm among EA is the Genetic Algorithm which was first introduced by [53]. The main mechanism behind GA is the ability to mimic biological evolution. GA is a population-based algorithm, meaning that it creates an initial set of random solutions called "Chromosomes" composed of genes as the elements of the solution; then, these solutions are subject to a "Fitness Function" defined by the user and, in our case is the performance of the SVR. To obtain the best solution for the optimization, GA used what are called genetic operators [54], that keep interacting with the solution until it converges toward the global optimum.
The SVR-GR model exhibits superior accuracy compared to the ANN model due to the utilization of a robust Genetic Algorithm, effectively reducing the error between the actual and predicted outputs. The optimized values of the SVR hyperparameters, obtained through the GR process, are shown in Table 3. Additionally, the parameters used in the optimization of the flow rate prediction through the use of SVR are outlined in Table 4 for further clarity. The R-score of the SVR-GA was 99.80%, as shown in Figure 10, with the root mean square error of 0.0149, and the mean absolute error of 0.0104 (see Figure 11).

Parameter
Value The size of the population 30 Maximum number of generations 8 Crossover probability 100 Mutation probability 20 The R-score of the SVR-GA was 99.80%, as shown in Figure 10, with the root mean square error of 0.0149, and the mean absolute error of 0.0104 (see Figure 11).  The R-score of the SVR-GA was 99.80%, as shown in Figure 10, with the root mean square error of 0.0149, and the mean absolute error of 0.0104 (see Figure 11).  Both developed data-driven models give powerful results with R scores higher than 98%. These models can be used to predict the flow rate for any new inputs with changes in the Fishbone well design parameters.

Models Comparison and Validation
In this section, we aim to perform a comprehensive comparison between the models discussed in Section 5. Our objective is to evaluate the performance and accuracy of each model, and to determine the most suitable model for our particular application.

Number of Branches
The results show that the flow rate increases with the increase in the number of branches. The trend and the results between the numerical simulation and the modified correlation are similar, with a relatively small difference compared to the modified analytical model (see Figure 12). The augmentation of the number of branches increases the exposed area of the reservoir to enhance production.

Number of Branches
The results show that the flow rate increases with the increase in the number of branches. The trend and the results between the numerical simulation and the modified correlation are similar, with a relatively small difference compared to the modified analytical model (see Figure 12). The augmentation of the number of branches increases the exposed area of the reservoir to enhance production.

Length of Branches
The results show that with the length of branches increasing, the flow rate increases in the numerical simulation, the modified correlation, and the analytical model, with a slight difference between them but nearly the same slope as shown in Figure 13. The reason for the production increase is the fact that a large portion of the reservoir has been exposed to the wellbores due to the length increase of the open hole branches, which aligns with the extended-reach drilling technology, with the idea to push the horizontal length to maximize reservoir drainage and minimize water production [55].

Length of Branches
The results show that with the length of branches increasing, the flow rate increases in the numerical simulation, the modified correlation, and the analytical model, with a slight difference between them but nearly the same slope as shown in Figure 13. The reason for the production increase is the fact that a large portion of the reservoir has been exposed to the wellbores due to the length increase of the open hole branches, which aligns with the extended-reach drilling technology, with the idea to push the horizontal length to maximize reservoir drainage and minimize water production [55].
The results show that the flow rate increases with the increase in the number of branches. The trend and the results between the numerical simulation and the modified correlation are similar, with a relatively small difference compared to the modified analytical model (see Figure 12). The augmentation of the number of branches increases the exposed area of the reservoir to enhance production.

Length of Branches
The results show that with the length of branches increasing, the flow rate increases in the numerical simulation, the modified correlation, and the analytical model, with a slight difference between them but nearly the same slope as shown in Figure 13. The reason for the production increase is the fact that a large portion of the reservoir has been exposed to the wellbores due to the length increase of the open hole branches, which aligns with the extended-reach drilling technology, with the idea to push the horizontal length to maximize reservoir drainage and minimize water production [55].

Permeability Anisotropy
The results found that permeability anisotropy affects productivity, as shown in Figure 14. When the anisotropy increases, the productivity increases. The numerical simulation and both developed analytical and empirical simulation showed a similar amount of production with the anisotropy variation. This means that intersecting the regions of higher permeability anisotropy will increase productivity. This finding supports the findings of [56], which studied the effect of permeability anisotropy on gas production, and as a result, reported that the higher horizontal permeability is more favorable for gas production. tion and both developed analytical and empirical simulation showed a similar amount of production with the anisotropy variation. This means that intersecting the regions of higher permeability anisotropy will increase productivity. This finding supports the findings of [56], which studied the effect of permeability anisotropy on gas production, and as a result, reported that the higher horizontal permeability is more favorable for gas production. Figure 14. The effect of permeability anisotropy on cumulative production.

Distance between Adjacent Branches
When the distance increases between the branches, the stimulated reservoir decreases, and the flow rate decreases slightly (see Figure 15). The numerical outputs show the same results as the developed correlation, with a constant difference compared to the developed analytical model. The three models confirm a similar trend of productivity decrease as a function of the distance. This validates the effect of the neighboring branches on each other based on the spacing between them, in the case where the distance is above the drainage radius. The results also demonstrate that the severity of the interference between the branches is very low with a very slight communication effect between them. The statistical comparison between the models as shown in Table 5 demonstrates their potential capabilities for productivity prediction of fishbone well design with some differences among them in terms of the coefficient of correlation and the average error. The developed analytical model presented low accuracy compared to the other models as Figure 14. The effect of permeability anisotropy on cumulative production.

Distance between Adjacent Branches
When the distance increases between the branches, the stimulated reservoir decreases, and the flow rate decreases slightly (see Figure 15). The numerical outputs show the same results as the developed correlation, with a constant difference compared to the developed analytical model. The three models confirm a similar trend of productivity decrease as a function of the distance. This validates the effect of the neighboring branches on each other based on the spacing between them, in the case where the distance is above the drainage radius. The results also demonstrate that the severity of the interference between the branches is very low with a very slight communication effect between them.
higher permeability anisotropy will increase productivity. This finding supports the find ings of [56], which studied the effect of permeability anisotropy on gas production, and as a result, reported that the higher horizontal permeability is more favorable for gas pro duction.

Distance between Adjacent Branches
When the distance increases between the branches, the stimulated reservoir de creases, and the flow rate decreases slightly (see Figure 15). The numerical outputs show the same results as the developed correlation, with a constant difference compared to the developed analytical model. The three models confirm a similar trend of productivity de crease as a function of the distance. This validates the effect of the neighboring branche on each other based on the spacing between them, in the case where the distance is above the drainage radius. The results also demonstrate that the severity of the interference be tween the branches is very low with a very slight communication effect between them. The statistical comparison between the models as shown in Table 5 demonstrate their potential capabilities for productivity prediction of fishbone well design with some differences among them in terms of the coefficient of correlation and the average error The developed analytical model presented low accuracy compared to the other models a Figure 15. The effect of the distance between adjacent branches on the cumulative production.
The statistical comparison between the models as shown in Table 5 demonstrates their potential capabilities for productivity prediction of fishbone well design with some differences among them in terms of the coefficient of correlation and the average error. The developed analytical model presented low accuracy compared to the other models as it is based on mathematical equations derived from Darcy's equation with many assumptions on the flow regime [57], compared to very sophisticated numerical methods to solve the complicated fundamental differential equation. The data-driven models had the advantages in terms of their high accuracy and low error because there were built in error terms, limited to the extent of the data collected; the same applied for the empirical model, which depends on the studied case and does not elucidate any underlying physics.

Conclusions
In this study, the authors developed three models (physical-based, empirical, and datadriven) to evaluate the gas well deliverability using Fishbone well drilling technology. The developed analytical model takes into consideration all Fishbone well design parameters, including the number, length of branches, distance between branches, and permeability anisotropy of the reservoir, offering superior features compared to existing models. It also has potential applications in unconventional reservoirs.
The comparison between the analytical and empirical models demonstrates their ability to accurately predict the gas flow rate well productivity by utilizing various parameters including the Fishbone well parameters, bottom hole pressures, and reservoir pressures. The results obtained from both models are consistent, indicating their effectiveness in this regard.
To improve the accuracy of analytical and empirical models, data-driven models were derived to help in effectively generating results with higher computational speed.
The three models can independently evaluate the well productivity and are useful for comparing results. Optimized with a Genetic Algorithm, the Support Vector Regression model showed the highest accuracy among the developed models. The applicability of these models depends on data availability, reservoir properties, and computational requirements.
Overall, the new analytical, empirical, and data-driven models for estimating the Fishbone well productivity are feasible and practical. The development and implementation of new productivity models for Fishbone well technology is imperative to advance the understanding of this cutting-edge drilling approach and its wide-ranging applications. This novel technology promises to revolutionize the way we recover hydrocarbons, elevating the level of precision and control in stimulation, while minimizing the associated hazards and costs. Furthermore, this innovative approach promotes a more sustainable and environmentally conscious method for enhanced recovery, making it an essential part of the industry's future.
Author Contributions: Each author contributed to the present paper. H.O. was responsible for preparing the methodology, analyzing the data, validating, drafting the paper, and drawing the conclusion. A.H. was responsible for the numerical simulation and data analyzation and validation. A.L. was responsible for reporting from the literature, analyzing the results and modifying the analytical model. A.C. was responsible for developing the data-driven models and reporting the results. V.R. reviewed the paper and discussed the results and methodology. M.M. provided the data and discussed the results. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author. The data are not publicly available due to confidentiality by the operator.