Predicting the Performance of Undeveloped Multi-Fractured Marcellus Gas Wells Using an Analytical Flow-Cell Model (FCM)

: The objective of the present study is to predict how changes in the fracture treatment design parameters will affect the production performance of new gas wells in a target zone of the Marcellus shale. A recently developed analytical ﬂow-cell model can estimate future production for new wells with different completion designs. The ﬂow-cell model predictions were benchmarked using historic data of 11 wells and 6 different completion designs. First, a type well was generated and used with the ﬂow-cell model to predict the performance of the later inﬁll wells—with variable completion designs—based off the performance of earlier wells. The ﬂow-cell model takes into account known hyperbolic forecast parameters ( q i , D i , and b -factor) and fracture parameters (height, half-length, and spacing) of a type well. Next, the ﬂow-cell model generates the hyperbolic decline parameters for an offset well based on the selected changes in the fracture treatment design parameters. Using a numerical simulator, the ﬂow-cell model was veriﬁed as an accurate modeling technique for forecasting the production performance of horizontal, multi-fractured, gas wells.


Introduction
This study uses detailed well data from a lease region in the Marcellus shale to benchmark history matching results of a recently developed analytical flow-cell model [1][2][3] and independent results from a numerical reservoir simulator (CMG-IMEX). The lease region in the Marcellus shale has 11 active gas production wells in a key part of the Middle Devonian organic-rich formation with ultra-low matrix permeability [4]. The Marcellus shale covers 95,000 square miles from New York State through Pennsylvania, West Virginia, and Ohio down to Kentucky [4]. Gas production from the Marcellus shale increased rapidly since 2008 when horizontal drilling and multistage hydraulic fracturing allowed economic recovery by increasing the reservoir contact and creating increased permeability [5]. The Marcellus shale exploitation has continued to grow and now is the biggest producer of natural gas within the United States at~23 Bcf/day out of the total production of~73 Bcf/day [6].
Throughout the years, operations within the Marcellus have changed significantly, driven by a variety of factors such as volatile commodity prices, varying well costs, and continual technology advancements. For example, the average lateral length was 4000 ft in 2011 [7] and by 2019 had more than doubled to 8500 ft, with some wells being drilled with lateral lengths over 15,000 ft [8]. The completion designs have also changed significantly. The average well in 2011 had a 350 ft stage spacing, and 1300 lbs of sand pumped per lateral foot [7]. Wells    Unfortunately, none of the traditional DCA methods can predict the performance of future wells with different completion designs, since it is an empirical fit of data that already exists. Until now, DCA methods can at best propose a type well that provides analogy performance, assuming the new well in the same formation will produce the same amount as the type well. In contrast, the flow-cell model proposed by Weijermars and coworkers [1][2][3] can predict future well performance using a type well, even when completion designs and well spacing are varied (see Section 2.2 and onward). Nelson et al. [12] used the three most popular decline curve equations (Arps, Duong, and Power Law Exponential) on field data from the Marcellus that were first history matched with a numerical simulator. With 30 years of production data generated by the numerical simulator, any of the three DCA methods turned out rather accurate ( Figure 2, R 2 ≥ 0.95). However, with less production history on a well, the various DCA methods tended to overpredict future production as they are not constrained by the physics of the reservoir. Nelson et al. [12] further found that the Arps equation gave the best fit for the long-term production behavior ( Figure 2). One must remember that the Power Law Exponential (PLE) method was aimed at mitigating some limitations of the single-segment Arps' equation by having a variable b-factor [13]. Duong [14] created an equation that is based on production dominated by fracture flow with negligible matrix contribution. The EURs using the Duong and PLE methods tend to be more conservative than those using Arps forecasts, especially in the later years of production [15]. Likewise, the Stretched Exponential Production Decline (SEPD), which was proposed to try model the performance of unconventional gas wells more accurately [16], did not evaluate favorably as compared to Arps [15].
Energies 2021, 14, x FOR PEER REVIEW 3 of 43 R 2 ≥ 0.95). However, with less production history on a well, the various DCA methods tended to overpredict future production as they are not constrained by the physics of the reservoir. Nelson et al. [12] further found that the Arps equation gave the best fit for the long-term production behavior ( Figure 2). One must remember that the Power Law Exponential (PLE) method was aimed at mitigating some limitations of the single-segment Arps' equation by having a variable b-factor [13]. Duong [14] created an equation that is based on production dominated by fracture flow with negligible matrix contribution. The EURs using the Duong and PLE methods tend to be more conservative than those using Arps forecasts, especially in the later years of production [15]. Likewise, the Stretched Exponential Production Decline (SEPD), which was proposed to try model the performance of unconventional gas wells more accurately [16], did not evaluate favorably as compared to Arps [15]. Typical multi-segmented decline curve for a shale well [5].

Figure 2.
Comparison of various decline curve predictions on a Marcellus well [12].
Unfortunately, none of the traditional DCA methods can predict the performance of future wells with different completion designs, since it is an empirical fit of data that already exists. Until now, DCA methods can at best propose a type well that provides analogy performance, assuming the new well in the same formation will produce the same amount as the type well. In contrast, the flow-cell model proposed by Weijermars and coworkers [1][2][3] can predict future well performance using a type well, even when completion designs and well spacing are varied (see Section 2.2 and onward).  [12].
Unfortunately, none of the traditional DCA methods can predict the performance of future wells with different completion designs, since it is an empirical fit of data that already exists. Until now, DCA methods can at best propose a type well that provides analogy performance, assuming the new well in the same formation will produce the Energies 2021, 14, 1734 4 of 42 same amount as the type well. In contrast, the flow-cell model proposed by Weijermars and co-workers [1][2][3] can predict future well performance using a type well, even when completion designs and well spacing are varied (see Section 2.2 and onward).

RTA Models for Marcellus Wells
Rate-Transient Analysis (RTA) also provides physics-based well performance solutions that incorporate different reservoir geometries with different reservoir flow regimes. While RTA models are not as quick as DCA models, they are generally fairly quick (<1 h per well) and sufficiently accurate, with the added benefit of maintaining a strong tie to the reservoir and physics compared to decline curve analysis. In unconventional reservoirs, RTA focuses on analyzing the transient linear flow regime with historical production data to estimate important well and reservoir parameters [5]. The transient linear flow solution for a horizontal, multi-fractured gas well is: 8T The left side of the equation is the normalized flowing pseudopressure for a gas well, made up of the flow rate (q), the initial pseudopressure (m(P i )), and bottom hole flowing pseudopressure (m(P wf )). The right side of the equation is made up of temperature (T), number of fractures (n f ), fracture height (H), fracture half-length (x f ), matrix permeability (k), porosity (φ), viscosity (µ g ), compressibility (c t ), and time (t). The The important unknowns in Equation (2) are the number of fractures, fracture height, fracture half-length, and matrix permeability while the rest of the variables can be estimated from logs, or cores. Unfortunately, there are four unknowns-but only one equationmeaning that the individual values cannot be solved. The combination of number of fractures, fracture height, fracture half-length, and fracture permeability are referred to as "Asqrtk" as shown in Equation (3). A simplified plot of normalized flowing pseudopressure versus the square root of time shows how to determine the slope, m, and the intercept, b' (Figure 3): Clearly, Asqrtk is a measure for productivity potential of a multi-fractured well-with large values indicating better drainage opportunity for the fluid in the reservoir-for a given viscosity and other fixed parameters for the reservoir and fracture properties.
Gunaydin et al. [5] developed an iterative workflow to estimate Asqrtk with the goals of honing in on the optimum fracture half-length, while varying the number of fractures. RTA was used to estimate the time of decline on a rate-constrained well. The model predicted the actual start of decline within 3 days accuracy, and the RTA-based production forecasts matched well on the production data after decline. In addition, on a different well, the EUR of the matched analytical model was within 3% of the numerical model. Their RTA case study was successful in predicting both the short term and long term production performance of the well.
Later studies used RTA to compare different completion designs in the Marcellus [17], using one pad with a new completion design and two pads with an older completion Energies 2021, 14, 1734 5 of 42 design. Looking at the production data, the new completion design had 20% higher production at 1 year. However, using RTA, the fracture contact area (Asqrtk) was 35% lower than the older completion design. Once the data was normalized for lateral length, bottom hole pressure profile, rock properties, and fluid properties, the analytical model for the new completion design only showed 8% more production after 1 year, while the EUR had reduced by 38%, due to the smaller fracture contact area. The study was successful in using analytical models to match historical production and estimate future production. Gunaydin et al. [5] developed an iterative workflow to estimate Asqrtk with the goals of honing in on the optimum fracture half-length, while varying the number of fractures. RTA was used to estimate the time of decline on a rate-constrained well. The model predicted the actual start of decline within 3 days accuracy, and the RTA-based production forecasts matched well on the production data after decline. In addition, on a different well, the EUR of the matched analytical model was within 3% of the numerical model. Their RTA case study was successful in predicting both the short term and long term production performance of the well.
Later studies used RTA to compare different completion designs in the Marcellus [17], using one pad with a new completion design and two pads with an older completion design. Looking at the production data, the new completion design had 20% higher production at 1 year. However, using RTA, the fracture contact area (Asqrtk) was 35% lower than the older completion design. Once the data was normalized for lateral length, bottom hole pressure profile, rock properties, and fluid properties, the analytical model for the new completion design only showed 8% more production after 1 year, while the EUR had reduced by 38%, due to the smaller fracture contact area. The study was successful in using analytical models to match historical production and estimate future production.
Trumbo et al. [18] presented a multivariate study of the completion design and its impact on well performance in the Utica (within the same basin as the Marcellus). Based on the unconventional transient linear flow analytical model, a predictive RTA tool was created using the Asqrtk per lateral length (Asqrtk/ft). Although the model satisfactory predicted a single variable (like Asqrtk/ft), the predictive model lacked a direct prediction of the resultant production profile necessary for economic analysis.

Recent Reservoir Simulations of Marcellus Wells
Filchock et al. [19] at the Marcellus Shale Energy and Environment Laboratory (MSEEL) tried to determine the impact of different completion parameters on reservoir performance. A numerical simulation based on 4 Marcellus wells in Northern West Virginia within the dry gas window found that increasing the fracture half-length by 100% only increased the EUR by 50%. MSEEL also analyzed those 4 wells in detail. Using various logs, there were over 1600 healed natural fractures found along a 6000 ft lateral, meaning there was a natural fracture on average every 4 ft [20]. The reservoir model of Filchock et al. [19] was able to match four years of historic production on two wells (unfortunately, Trumbo et al. [18] presented a multivariate study of the completion design and its impact on well performance in the Utica (within the same basin as the Marcellus). Based on the unconventional transient linear flow analytical model, a predictive RTA tool was created using the Asqrtk per lateral length (Asqrtk/ft). Although the model satisfactory predicted a single variable (like Asqrtk/ft), the predictive model lacked a direct prediction of the resultant production profile necessary for economic analysis.

Recent Reservoir Simulations of Marcellus Wells
Filchock et al. [19] at the Marcellus Shale Energy and Environment Laboratory (MSEEL) tried to determine the impact of different completion parameters on reservoir performance. A numerical simulation based on 4 Marcellus wells in Northern West Virginia within the dry gas window found that increasing the fracture half-length by 100% only increased the EUR by 50%. MSEEL also analyzed those 4 wells in detail. Using various logs, there were over 1600 healed natural fractures found along a 6000 ft lateral, meaning there was a natural fracture on average every 4 ft [20]. The reservoir model of Filchock et al. [19] was able to match four years of historic production on two wells (unfortunately, without mentioning if natural fractures were included in the model).
Several studies have suggested the importance of stress-dependent permeability and conductivity within the Marcellus [9,21]. The studies showed an order of magnitude reduction in both the matrix permeability and the fracture conductivity, due to net stress increases when the reservoir pressure is depleted. El Sgher et al. [22] used a reservoir simulator on two Marcellus wells to investigate the impact of the individual hydraulic fractures. The principal goal was to understand geomechanical effects of the hydraulic fractures and non-uniform formation properties on the production and gas recovery with horizontal, multi-fracture wells in the Marcellus Shale. The study of El Sgher et al. [22] came about since multiple sources had suggested that only 50-60% of the perforation clusters were properly stimulated and producing at meaningful rates in other shale resource plays.
El Sgher et al. [22] history matched the production profile with four plausible solutions by varying the fracture half-length and the hydraulic fracture model. The first numerically simulated well had no natural fractures, the second had natural fractures every 20 ft and the last two had a different amount of "missing" hydraulic fractures (50% reduction in halflength of the central-stress-shadowed-perforation clusters within a pumped stage with high ISIP and high clay content). All of the models had the stress-dependent permeability and conductivity. The matches of the first two models overpredicted the first four months of production. The fracture half-lengths did not have to be varied to achieve a good fit in the first model, while the authors claimed that the hydraulic fractures are propped and draining 400 ft of height (even though the Marcellus target zone is only 100 ft thick). The fracture half-lengths had to significantly reduce (15%) when the reservoir was modeled with natural fractures. The naturally fractured model had a 25% increase in EUR. The authors did not show actual simulation results for the last two models with stages missing, but nonetheless concluded that the fracture stages do not equally contribute to production.
A major conclusion of El Sgher et al. [22] was that not all hydraulic fractures contribute uniformly to production and it is important to determine the distribution of the natural fractures along the wellbore length. In addition, their fracture model predicted fracture half-lengths and fracture conductivity in a reservoir without natural fractures. Some of the inputs needed for the reservoir numerical simulations in the present study (see Section 4) were based off the study by El Sgher et al. [22].
Chin et al. [23] presented a practical workflow to predict child well performance utilizing a numerical simulator that included pressure depletion on initial reservoir pressure and incorporated geomechanical effects on cluster efficiency. When initially trying to match the production of the child well, no combination of fracture half-length, conductivity, and matrix perm matched early time production data. Reducing the initial pressure brought down the early time production and created a better match. There should not have been any pressure depletion from fluid flow within the ultra-low permeability matrix. This led the authors to conclude that the fractures networks were connected between the parent and child wells. Due to the connected fracture networks, the cluster efficiency was reduced from 50% to 40% based on the assumption that new fractures would not propagate once a stage became connected to an existing fracture network. The reduction of the cluster efficiency was effective at getting a better match on production data. On the other hand, Weijermars et al. [24] encountered similar problems in pressure matching well data due to unconstrained hydraulic fracture height, which can be solved when the fracture terminates against the upper and lower boundaries of the producing target zone.
Deasy et al. [25] presented a method of testing and confirming the production results of varying completion designs within the liquid rich portion of the Marcellus starting with a numerical simulator. The first step was using a numerical simulator to match the production of a certain completion design (300 ft stage spacing). Next, fracture stage spacing was varied to find the optimal rate of return on the investment. The new completion design (150 ft stage spacing) was tested in a field trial. An important part of the field testing was control for geologic properties, interior versus exterior, lateral length, well spacing, similar line pressures, vintage, and completion properties per stage (lb/stg and bbl/stg). The field test proved that that the numerical simulator gave reliable predictions, with the new wells showing a 60% increase in initial production and a 20% increase in the EUR.
However, for practical purposes, the operator needed a quick way to create a forecast or type well in areas outside of the pilot test [25]. From the results of the numerical simulator, a so-called (production) uplift ratio was abstracted from the comparison of the two (early and modified) completion designs ( Figure 4). The uplift ratio sacrificed some accuracy by not being fully tied to the numerical simulator, but a blind test of the uplift ratio showed that the forecasted volumes were accurate within 5% of the actual volumes produced. A full implementation-by the operator-of the new completion design, followed.
Instead of resorting to the uplift ratio for forecasting future well results, the flow-cell model [1] keeps the forecast tied to a type well, which could also be paired with a numerical reservoir simulator for benchmarking, as shown in the present study. (early and modified) completion designs ( Figure 4). The uplift ratio sacrificed some accuracy by not being fully tied to the numerical simulator, but a blind test of the uplift ratio showed that the forecasted volumes were accurate within 5% of the actual volumes produced. A full implementation-by the operator-of the new completion design, followed.
Instead of resorting to the uplift ratio for forecasting future well results, the flow-cell model [1] keeps the forecast tied to a type well, which could also be paired with a numerical reservoir simulator for benchmarking, as shown in the present study.

Flow-Cell Based Models-A New Approach
Weijermars and Khanal [1] presented a physics-based flow-cell model for estimating the production of new, undeveloped wells. The complex analysis method (CAM) is an analytical solution that has high resolution, while allowing for quick calculations of streamlines in a flow cell. The flow-cell model is built upon describing the flow paths of a fluid particles in the reservoir space between two adjacent hydraulic fractures confined vertically with no flow boundaries above and below [1]. The flow cell is defined by D (fracture spacing), H (fracture height), and L (fracture length) as shown in Figure 5a. Multiple flow-cells can be combined to represent a horizontal, multi-fractured well where each flow cell-in the basic model-is assumed to be the same (Figure 5b).

Flow-Cell Based Models-A New Approach
Weijermars and Khanal [1] presented a physics-based flow-cell model for estimating the production of new, undeveloped wells. The complex analysis method (CAM) is an analytical solution that has high resolution, while allowing for quick calculations of streamlines in a flow cell. The flow-cell model is built upon describing the flow paths of a fluid particles in the reservoir space between two adjacent hydraulic fractures confined vertically with no flow boundaries above and below [1]. The flow cell is defined by D (fracture spacing), H (fracture height), and L (fracture length) as shown in Figure 5a. Multiple flow-cells can be combined to represent a horizontal, multi-fractured well where each flow cell-in the basic model-is assumed to be the same (Figure 5b). duced. A full implementation-by the operator-of the new completion design, followed.
Instead of resorting to the uplift ratio for forecasting future well results, the flow-cell model [1] keeps the forecast tied to a type well, which could also be paired with a numerical reservoir simulator for benchmarking, as shown in the present study.

Flow-Cell Based Models-A New Approach
Weijermars and Khanal [1] presented a physics-based flow-cell model for estimating the production of new, undeveloped wells. The complex analysis method (CAM) is an analytical solution that has high resolution, while allowing for quick calculations of streamlines in a flow cell. The flow-cell model is built upon describing the flow paths of a fluid particles in the reservoir space between two adjacent hydraulic fractures confined vertically with no flow boundaries above and below [1]. The flow cell is defined by D (fracture spacing), H (fracture height), and L (fracture length) as shown in Figure 5a. Multiple flow-cells can be combined to represent a horizontal, multi-fractured well where each flow cell-in the basic model-is assumed to be the same (Figure 5b). The particle flow paths along with the particle velocities are used to calculate the flow rate and pressures in a flow cell over time. The flow rate and pressure of a single flow cell is then combined with all of the other flow cells to represent the whole well. The pressure depletion rate is then scaled based off a type well. With the assumption that bottom hole pressure is the same between the type well and a new well, and allowing for design changes in the dimension of the flow cells for the new wells (fracture spacing, fracture height, fracture half-length, and total well-length), it is possible to determine the flow rate of a new well based on the flow performance of the type well. Besides assuming that both the new well and the type well have the same BHP, the flow-cell model makes the following simplifying assumptions:

•
There are the same reservoir and fracture properties between all flow cells as they are all considered identical.

•
The proppant concentration over the considered fracture half-length is high enough to ensure near infinite conductivity meaning that pressure drop within the fractures is negligible.

•
There is no complex fracture geometry. All flow cells are bound by planar fractures, transverse to the wellbore.

•
The impact of flow from the unconfined matrix to the first fracture and last fracture of the well is neglected by assuming the same boundary conditions for all flow cells.

•
There are similar production operations between the type well and the new well and no surface flow constrictions.

•
The production impact outside of the flow-cell is not neglected, but is assumed to give flow contributions to the fracture tips in both the type well and the new well ( Figure 5c).

•
There may be an overall change in hydrocarbon recovery percentage over the life of the new well, but only when the flow-cell dimensions increase as compared to the type well. However, the effect is mostly due to accelerated early production when fracture spacing is reduced. • Gains in EUR would come from increasing the flow-cell dimensions, which is possible by (1) increasing the fracture height, (2) increasing the fracture half-length, or (3) increasing the well length.

•
The well performance predicted by the flow-cell model is captured here by a singlesegment hyperbolic decline curve.
When the limiting conditions are fulfilled, the flow-cell model can be used. The CAM flow model starts with a known Arps production forecast (whether it is a single well or type well) and then scales the flow rates based on the completion design of the new well. The particle flow paths along with the particle velocities are used to calculate the flow rate and pressures in a flow cell over time. The flow rate and pressure of a single flow cell is then combined with all of the other flow cells to represent the whole well. The pressure depletion rate is then scaled based off a type well. With the assumption that bottom hole pressure is the same between the type well and a new well, and allowing for design changes in the dimension of the flow cells for the new wells (fracture spacing, fracture height, fracture half-length, and total well-length), it is possible to determine the flow rate of a new well based on the flow performance of the type well. Besides assuming that both the new well and the type well have the same BHP, the flow-cell model makes the following simplifying assumptions:

•
There are the same reservoir and fracture properties between all flow cells as they are all considered identical.

•
The proppant concentration over the considered fracture half-length is high enough to ensure near infinite conductivity meaning that pressure drop within the fractures is negligible.

•
There is no complex fracture geometry. All flow cells are bound by planar fractures, transverse to the wellbore.

•
The impact of flow from the unconfined matrix to the first fracture and last fracture of the well is neglected by assuming the same boundary conditions for all flow cells.

•
There are similar production operations between the type well and the new well and no surface flow constrictions.

•
The production impact outside of the flow-cell is not neglected, but is assumed to give flow contributions to the fracture tips in both the type well and the new well ( Figure 5c).

•
There may be an overall change in hydrocarbon recovery percentage over the life of the new well, but only when the flow-cell dimensions increase as compared to the type well. However, the effect is mostly due to accelerated early production when fracture spacing is reduced. • Gains in EUR would come from increasing the flow-cell dimensions, which is possible by (1) increasing the fracture height, (2) increasing the fracture half-length, or (3) increasing the well length.

•
The well performance predicted by the flow-cell model is captured here by a singlesegment hyperbolic decline curve.
When the limiting conditions are fulfilled, the flow-cell model can be used. The CAM flow model starts with a known Arps production forecast (whether it is a single well or type well) and then scales the flow rates based on the completion design of the new well. The model uses fracture dimensions (fracture half-length, number of fractures, and fracture spacing) to scale the flow rate. The practical value of the model is that it provides the Energies 2021, 14, 1734 9 of 42 hyperbolic forecast with the Arps parameters (q i , D i , and b-factor) for the new well using solely changes in the well design parameters.
The flow-cell model is made up of 5 equations to estimate the Arps production forecast inputs for a new well. The first equation of the flow-cell model determines the initial flow rate of the new well (q i_NewWell ) and is shown in Equation (4). The initial flow rate of the new well is calculated based on the initial flow rate of the type well (q i_Typewell ), the fracture spacing of the new well (D N ), the fracture spacing of the type well (D T ), the number of fractures of the new well (N N ), the number of fractures of the type well (N T ), and the ratio of the drained volume of new well and type well (V R ). The initial flow rate of the new well is sensitive to the change in number of fractures and fracture spacing between the new well and the type well.
The initial decline rate of the new well (D i_N ) is shown in Equation (5). The initial decline rate of the new well is calculated based on the initial decline rate of the type well (D i_T ), the fracture spacing of the new well (D N ), the fracture spacing of the type well (D T ), the number of fractures of the new well (N N ), the number of fractures of the type well (N T ), and the ratio of the drained volume of new well and type well (V R ). The initial decline rate of the new well is very sensitive to the change in fracture spacing between the new well and the type well.
The b-factor for the new well (b N ) is determined using Equation (6). The b-factor of the new well is calculated based on the b-factor of the type well (b T ), the initial decline rate of the new well (D i_N ), the initial decline rate of the type well (D i_T ), the fracture spacing of the new well (D N ), the fracture spacing of the type well (D T ), the number of fractures of the new well (N N ), the number of fractures of the type well (N T ), and the ratio of the drained volume of new well and type well (V R ).
The volume ratio (V R ) is determined using Equation (7). The volume ratio accounts for any differences in the drainage areas between the new well and the type well. The volume ratio is calculated based on the fracture spacing of the new well (D N ), the fracture spacing of the type well (D T ), the number of fractures of the new well (N N ), the number of fractures of the type well (N T ), the length of the fractures of the new well (L N ), the length of the fractures of the type well (L T ), the fracture height of the new well (H N ), and the fracture height of the type well (H T ).
With the initial flow rate, initial decline rate, and b-factor of the new well determined, the flow rate of the new well at any point in time (q(t)) can be determined using the Arps equation as shown in Equation (8).
For more details on the flow-cell model derivation please see the original papers by Weijermars and coauthors [1][2][3].
Utilizing the flow-cell model, Weijermars and Khanal [1] were able to effectively and efficiently predict future well performance for various fracture spacing, based on the type well performance. As shown in Figure 6, the oil rate and cumulative oil production were calculated over time for a number undeveloped wells in the oil window of the Eagle Ford Formation (East Texas), considering various fracture spacing. A similar analysis will be performed in the present paper for undeveloped gas wells in the Marcellus Formation. After first benchmarking and verifying the flow-cell model with a numerical simulator, the model is next applied to field data of 11 legacy wells (see Sections 3 and 4).
For more details on the flow-cell model derivation please see the original papers by Weijermars and coauthors [1][2][3].
Utilizing the flow-cell model, Weijermars and Khanal [1] were able to effectively and efficiently predict future well performance for various fracture spacing, based on the type well performance. As shown in Figure 6, the oil rate and cumulative oil production were calculated over time for a number undeveloped wells in the oil window of the Eagle Ford Formation (East Texas), considering various fracture spacing. A similar analysis will be performed in the present paper for undeveloped gas wells in the Marcellus Formation. After first benchmarking and verifying the flow-cell model with a numerical simulator, the model is next applied to field data of 11 legacy wells (see Sections 3 and 4).

Case study Well Data and Basic Production Analysis
An operator provided historic production data for 11 wells within the dry gas portion of the Marcellus. The well spacing, and completion timing of the 11 wells are explained, together with the completion designs pumped as well as the estimated fracture properties. The fracture dimensions are important inputs into the flow-cell model. The reservoir flow regimes are identified and the production history presented, as the flow rates and Arps decline parameters are key inputs for the flow-cell model.

Location of Case Study Wells
The focus of this study is on the Appalachian Basin specifically the dry gas window of the Marcellus (Middle Devonian formation) in Pennsylvania. A subset of 11 wells drilled from 4 well pads within the dry gas window will be studied in detail. All wells were landed in the lower Marcellus target zone, which occurs sandwiched between two limestone formations, i.e., the Cherry Valley (top) and the Onondoga (bottom). Both formations are assumed to act as fracture barriers, limiting the height of fracture growth into the Upper Marcellus [23]. A gun-barrel view of the wells with their horizontal spacing and the year of completion is shown in Figure 7. To protect the privacy of the company that graciously provided data for the study, the paper does not go into any further detail of the location of the wells.
The full schedule of well development is given in Table 1

Case study Well Data and Basic Production Analysis
An operator provided historic production data for 11 wells within the dry gas portion of the Marcellus. The well spacing, and completion timing of the 11 wells are explained, together with the completion designs pumped as well as the estimated fracture properties. The fracture dimensions are important inputs into the flow-cell model. The reservoir flow regimes are identified and the production history presented, as the flow rates and Arps decline parameters are key inputs for the flow-cell model.

Location of Case Study Wells
The focus of this study is on the Appalachian Basin specifically the dry gas window of the Marcellus (Middle Devonian formation) in Pennsylvania. A subset of 11 wells drilled from 4 well pads within the dry gas window will be studied in detail. All wells were landed in the lower Marcellus target zone, which occurs sandwiched between two limestone formations, i.e., the Cherry Valley (top) and the Onondoga (bottom). Both formations are assumed to act as fracture barriers, limiting the height of fracture growth into the Upper Marcellus [23]. A gun-barrel view of the wells with their horizontal spacing and the year of completion is shown in Figure 7. To protect the privacy of the company that graciously provided data for the study, the paper does not go into any further detail of the location of the wells.   N7  2014  2014  2015  S1  2015  2015  2015  N1  2014  2015  2015  N9  2015  2015  2015  S2  2015  2015  2015  N8  2015  2015  2015  N4  2015  2016  2016  N5  2015  2018  2018  N6  2015  2018  2018  N2  2015  2019  2019  N3  2015  2019  2019 The 11 wells of this study were all drilled in approximate N-S direction and are all parallel (Figure 7). The majority of the wells have a spacing of about 1300 ft, only the easternmost wells are spaced at 2500 ft, leaving possible infill space for another, future well in between. All the current wells target the same interval in the Lower Marcellus, which is around 100 ft thick (Figure 7).

Completion Parameters
The completion design for each of the wells varies, depending on the year of completion. Table 2 has more general completion information, while Table 3 has the information per stage, per cluster and per foot of completed lateral. Over time the number of perfora-  The full schedule of well development is given in Table 1. One well was drilled and completed in 2014 while all other 10 wells were drilled in 2015. Due to the gas price crash during 2015, not all of the wells were completed right after being drilled. Only 5 of the wells were completed in 2015, 1 in 2016, 2 more in 2018, and the last 2 in 2019. The 11 wells of this study were all drilled in approximate N-S direction and are all parallel ( Figure 7). The majority of the wells have a spacing of about 1300 ft, only the easternmost wells are spaced at 2500 ft, leaving possible infill space for another, future well in between. All the current wells target the same interval in the Lower Marcellus, which is around 100 ft thick (Figure 7).

Completion Parameters
The completion design for each of the wells varies, depending on the year of completion. Table 2 has more general completion information, while Table 3 has the information per stage, per cluster and per foot of completed lateral. Over time the number of perforation clusters has increased while the spacing between each cluster has significantly decreased (100 ft to 40 ft). The sand per foot of lateral has generally been around 2000 lb/ft while the amount of sand per cluster has decreased significantly as more clusters were added without adding more sand. The amount of fracturing fluid pumped per foot has also increased in the more recent designs, but has decreased slightly on a fluid per cluster base (Table 3).

Fracture Half-Lengths
IHS Harmony was used to evaluate the drainage area of each of the wells to understand the impact of changing the completion designs. Asqrtk, as mentioned in Section 2.1.3, is the amount of area that the fractures contact multiplied by the square root of permeability, which is a proxy for how effective the reservoir was stimulated by the fracturing. When a well is within linear flow, the boundaries of the reservoir are not known, but Asqrtk can still be determined. The superposition time and variable pressure models were used to obtain Asqrtk from the slope of a tangent line on the normalized pressure vs. square root time of plot (e.g., Figure 3).
Next, the fracture half-length was estimated using the values determined from the superposition time model. The fracture half-length, fracture conductivity, and permeability of the SRV were varied to minimize the error on the history match. Table 4 contains the resultant Asqrtk and the inferred fracture half-lengths for the case study wells. As the number of perforation clusters increased, the Asqrtk increased while the fracture half-length decreased. The trend becomes even more obvious when the completions are binned to modern completion designs (2016+ completion year) against the early completion designs (2014-2015) ( Table 5). The fracture half-length decreased by 17%, but the Asqrtk per foot of lateral increased by 27% between the early and modern completion designs.

Flow Regimes
Identifying flow regimes is important for obtaining improved history matches with DCA curves and may lower the uncertainty in the production forecasts and the estimated ultimate recovery (EUR) [11]. The Fetkovich type curve was used to identify flow regimes. All of the wells studied are multi-fractured wells, meaning they are initially expected to be in the linear flow regime which is a 1 2 slope on the Fetkovich type curve. If there is boundary-dominated flow (BDF), a unit slope or steeper appears on the type curve.
As seen in Figure 8, the newest wells, N2 and N3, exhibit a strong linear flow trend, but without meeting any boundaries since the slopes are 1/2 for the entire production history. Both the N2 and N3 wells had produced for around 1.5 years at the time of our analysis. Even the oldest wells (N1 and S1, Figure 8) after producing for almost 4.5 years did not yet exhibit BDF. In conclusion, none of the wells in our case study show-even after several years of production-any obvious sign of BDF.

Flow Regimes
Identifying flow regimes is important for obtaining improved history matches with DCA curves and may lower the uncertainty in the production forecasts and the estimated ultimate recovery (EUR) [11]. The Fetkovich type curve was used to identify flow regimes. All of the wells studied are multi-fractured wells, meaning they are initially expected to be in the linear flow regime which is a ½ slope on the Fetkovich type curve. If there is boundary-dominated flow (BDF), a unit slope or steeper appears on the type curve.
As seen in Figure 8, the newest wells, N2 and N3, exhibit a strong linear flow trend, but without meeting any boundaries since the slopes are 1/2 for the entire production history. Both the N2 and N3 wells had produced for around 1.5 years at the time of our analysis. Even the oldest wells (N1 and S1, Figure 8) after producing for almost 4.5 years did not yet exhibit BDF. In conclusion, none of the wells in our case study show-even after several years of production-any obvious sign of BDF.

Production History
To accurately predict future production, the decline curve analysis should be tied to the understanding of the production method and the production data. The case study

Production History
To accurately predict future production, the decline curve analysis should be tied to the understanding of the production method and the production data. The case study wells (Table 6) are operated with two different initial production methods: (1) All wells that started producing before 2018 had tubing installed right after completions; (2) The more recent wells had no tubing installed, and instead produce out of just the casing. The change in initial production method reduced the friction within the wellbore significantly, allowing for higher initial flow rates and drawdown in the newer wells. Figure 9 shows the production histories for each of the case study wells (colored by the year of the completion). Changes in the more recent completion designs (2016-2019) have resulted in increased early time production.  In addition, some of the parent wells experienced fracture hits. Any time a child well was completed next to a parent, the parent well was shut-in and had a few months without production, after which pressure build-up was needed for the parent well to unload the extra water.

Flow-Cell Model Workflow and Results
For the flow-cell model to be beneficial and used by reservoir engineers, it must save time and be easier to use as compared to a numerical simulator; however, it must also predict correctly the performance. This paper will now proceed to verify-using a numerical simulator-the accuracy of the flow-cell model as a proxy reservoir model for a dry gas reservoir. The flow-cell model has previously been verified in a black oil well in the Eagle Ford Formation [1][2][3]. However, the flow-cell model has never before been tested and verified for a dry gas well. The present study uses a CMG history match to benchmark the flow-cell model prediction and confirms that the flow-cell model is a valid tool for use in dry gas wells (see later). Subsequent to the model verification, the flow-cell model will be used to predict the performance of wells with different completion designs and then compare those predictions to the actual performance of the wells. The flow-cell model, thus first verified, was next used in a case study of actual well performance, with different completion designs and well spacing.
The adopted work flow is as follows: 1.
Generate single-segment Arps decline curves of the 11 wells using regression and engineering judgment.

2.
Compare and contrast the Arp decline curve inputs (q i , D i , and b-factor) for 5-year and 50-year cumulative production.

3.
Use the horizontal, multi-fracture analytical model to history match the pressures and rates to estimate the fracture half-lengths of each well.

4.
Use the flow-cell model to predict neighboring well performance based on the actual well completion parameters (fracture height and number of clusters) and the fracture half-length from the analytical model.

5.
Compare the flow-cell model prediction to actual production of the more recent wells. Compare production rate vs. time and cumulative production vs. time in a 5-year and 50-year timeframe.

Verifying the Reservoir Model
A commercial dry gas reservoir simulator (CMG-IMEX) was used to history match a 30-year gas production forecast for a type well created from five of the case study wells completed in 2015 ( Table 1). The starting point for the numerical simulation was to select model properties from prior literature. The model properties from Chin et al. [23] were preferentially used as that study also focused on the northeastern Pennsylvania portion of the Marcellus. Table 7 shows the model properties held constant for all of the numerical model realizations. The lateral spacing of the wells was 1300 ft (Figure 7; we ignore the delay in wells' early production). The width of the reservoir model for the type well was enlarged to 3000 ft (to allow for space beyond the fracture half-length and drained volume). The length of the reservoir model was also increased beyond the well length to avoid boundary effects. The thickness of the reservoir was not increased beyond the pay zone height of 107 ft, because there are fracture barriers (Cherry Valley and Onondoga Limestone) on either side of the Marcellus [23]. The type well was designed with a constant cluster spacing of 100 ft, which was kept constant in the history matched numerical model. With the well length of 6450 ft, the type well has 65 transverse hydraulic fractures. All of the fractures were assumed to be 100% open and contributing to the well productivity. The type well is assumed to have a constant bottom hole flowing pressure of 600 psi. Unknown reservoir and fracture properties were varied to find the optimal match of daily gas production and the cumulative gas production. The optimal match was subsequently used to verify the accuracy of the well rate forecast with the flow-cell model for various fracture spacing.
Some of the actual reservoir and fracture properties remain poorly constrained or unknown, like fracture half-length, fracture height, fracture width, fracture permeability, matrix permeability, matrix porosity, and initial water saturation. Those properties were allowed to vary throughout the history match to allow the numerical reservoir model to try to match the production forecast for the type well ( Figure 10). The properties of the initial guess and best match type well case are shown in Table 8. The type well forecast was used for all 30 years instead of using only the actual production data for the first four years. The initial rate was 5065 Mcf/day with an initial nominal decline rate of 1.009 and a b-factor of 1.13. . Daily gas production rate (blue) and cumulative gas production (orange) over 30 years for a type well with 0 months of historical data and 30 years of forecast, which was used as the basis of the history match for the numerical simulator. Figure 10. Daily gas production rate (blue) and cumulative gas production (orange) over 30 years for a type well with 0 months of historical data and 30 years of forecast, which was used as the basis of the history match for the numerical simulator. The sensitivity analysis and history match within CMOST were used to find the numerical model to best match the type well daily production and cumulative gas production. The sensitivity analysis used a random seed to generate a parameter from an input distribution. All of the parameters were given uniform distributions with large but practical bounds. The sensitivity analysis used randomized variables to history match and minimize the mismatch between the type well forecast and the numerical model. The default solver within CMG (Design Exploration Controlled Evolution optimization algorithm) was used. The solver used regression to try to predict the numerical model with the error minimized and then produced the forecasts for that model. The history match took over 200 different scenarios to find an optimal match. Some of the sensitivity analysis and history matching scenarios are shown in Figure 11. The optimal match of the type well forecast from the numerical model is shown in Figure 12. The optimal match had an average absolute relative error of 2.3% and 15.2% over 30 years for cumulative gas production and daily gas production, respectively. The sensitivity analysis and history match within CMOST were used to find the numerical model to best match the type well daily production and cumulative gas production. The sensitivity analysis used a random seed to generate a parameter from an input distribution. All of the parameters were given uniform distributions with large but practical bounds. The sensitivity analysis used randomized variables to history match and minimize the mismatch between the type well forecast and the numerical model. The default solver within CMG (Design Exploration Controlled Evolution optimization algorithm) was used. The solver used regression to try to predict the numerical model with the error minimized and then produced the forecasts for that model. The history match took over 200 different scenarios to find an optimal match. Some of the sensitivity analysis and history matching scenarios are shown in Figure 11. The optimal match of the type well forecast from the numerical model is shown in Figure 12. The optimal match had an average absolute relative error of 2.3% and 15.2% over 30 years for cumulative gas production and daily gas production, respectively.
(a) (b) Figure 11. Various sensitivity analysis and history matching scenarios: (a) daily gas production rate vs. time; (b) cumulative gas production vs. time. Figure 11. Various sensitivity analysis and history matching scenarios: (a) daily gas production rate vs. time; (b) cumulative gas production vs. time.
(a) (b) Figure 11. Various sensitivity analysis and history matching scenarios: (a) daily gas production rate vs. time; (b) cumulative gas production vs. time.
(a) (b) Figure 12. The optimal match of the numerical simulation and the type well production (Arps DCA): (a) daily gas production rate vs. time; (b) cumulative gas production vs. time. Figure 12. The optimal match of the numerical simulation and the type well production (Arps DCA): (a) daily gas production rate vs. time; (b) cumulative gas production vs. time.

Flow-Cell Model Benchmarked against Numerical Simulator
The numerical simulator was run with 3 different fracture spacing (50, 70, and 150 ft) from the optimal match of the type well (100 ft). The fracture spacing was the only variable changed from the optimal match model properties. It was assumed that all of the same fracture properties will be created just at a different spacing interval. The daily gas production and cumulative gas production for the numerical models with different fracture spacing are shown in Figure 13.

Flow-Cell Model Benchmarked against Numerical Simulator
The numerical simulator was run with 3 different fracture spacing (50, 70, and 150 ft) from the optimal match of the type well (100 ft). The fracture spacing was the only variable changed from the optimal match model properties. It was assumed that all of the same fracture properties will be created just at a different spacing interval. The daily gas production and cumulative gas production for the numerical models with different fracture spacing are shown in Figure 13.
The flow-cell model was used with the type well DCA (100 ft fracture spacing) to predict the performance for the cases of different fracture spacing (50, 70, and 150 ft). The fracture spacing was the only variable changed in the flow-cell model. The daily gas production and cumulative gas production for the initial 100 ft spacing and three different fracture spacing are shown in Figure 14a and b, respectively.
(a) (b) Figure 13. The optimal match of the numerical simulation with varying fracture spacing and the resultant (a) daily gas production rate and (b) cumulative gas production over 30 years. Figure 13. The optimal match of the numerical simulation with varying fracture spacing and the resultant (a) daily gas production rate and (b) cumulative gas production over 30 years.
The flow-cell model was used with the type well DCA (100 ft fracture spacing) to predict the performance for the cases of different fracture spacing (50, 70, and 150 ft). The fracture spacing was the only variable changed in the flow-cell model. The daily gas production and cumulative gas production for the initial 100 ft spacing and three different fracture spacing are shown in Figure 14a,b, respectively.
(a) (b) Figure 13. The optimal match of the numerical simulation with varying fracture spacing and the resultant (a) daily gas production rate and (b) cumulative gas production over 30 years.
(a) (b) Figure 14. The flow-cell model using the type well DCA with varying fracture spacing and the resultant (a) daily gas production rate vs. time and (b) cumulative gas production vs. time. The flow-cell model and the numerical simulation were based off of the same type well production forecast. Both models were run with the same varying fracture spacing. The cumulative gas production curves for the 4 different fracture spacing for both models are shown in Figure 15. The cumulative gas production for both models vary over time, and have 30-year EURs with similar values ( Table 9). The close match of the EURs from the numerical simulation results with those from the flow-cell model verifies that the flow-cell model can be confidently used to estimate the production for new wells when fracture treatment design parameters (such as fracture spacing) are changed. The flow-cell model and the numerical simulation were based off of the same type well production forecast. Both models were run with the same varying fracture spacing. The cumulative gas production curves for the 4 different fracture spacing for both models are shown in Figure 15. The cumulative gas production for both models vary over time, and have 30-year EURs with similar values ( Table 9). The close match of the EURs from the numerical simulation results with those from the flow-cell model verifies that the flowcell model can be confidently used to estimate the production for new wells when fracture treatment design parameters (such as fracture spacing) are changed. All of the cases had low error on the cumulative gas production at 30 years (under 5% mismatch in each scenario) as shown in Table 9. In both the numerical simulation and the flow-cell model, there is cumulative gas production benefit from decreasing fracture spacing as there is a ~10% increase of the 30-year EUR when the fracture spacing tightens  All of the cases had low error on the cumulative gas production at 30 years (under 5% mismatch in each scenario) as shown in Table 9. In both the numerical simulation and the flow-cell model, there is cumulative gas production benefit from decreasing fracture spacing as there is a~10% increase of the 30-year EUR when the fracture spacing tightens from 100 to 50 ft. Most of the benefit of decreasing the fracture spacing comes from the early production acceleration. The 50 ft fracture spacing model produced 35% more gas than the 100 ft fracture spacing model at 2,000 days. Table 10 compares the 30-year cumulative production of the new wells with different fracture spacing to the type well of 100 ft fracture spacing. The small difference between the numerical simulation and the flow-cell model shows that the flow-cell model is an accurate reservoir model.

Mismatch of Field Case and Verified Flow-Cell Forecasts
Using the flow-cell model [Equations (4)-(8) from Section 2.2], the initial flow rate, initial decline, and the b-factor were calculated for the six wells with different completion designs based on the type well, as shown in Table 11. The normalized initial flow rates from the flow-cell model were around 2 to 3 times larger than the least-squared regression initial flow rate. The initial nominal decline rates varied more and were around 1.5 to 6 times larger than the least-squared regression initial declines. The b-factor changed the least of the Arps decline parameters between the two methods with some wells seeing b-factor increases up to 25% and other wells see decreases up to 7%. The normalized 5-year and 50-year recoveries for the least-squared DCA regression and the flow-cell model are compared in Table 12. There is a significant amount of variance between the two techniques. Of the 6 wells, the flow-cell model largely overpredicted (greater than 5%) the 5-year recovery on 4 wells, slightly underpredicted (less than −5%) on 1 well, and correctly predicted (mismatch less than absolute value of 5%) on 1 well. The flow-cell model correctly predicted the 50-year recovery on 2 wells, overpredicted on 2 wells, and underpredicted on the remaining two wells. The well (N6) that was correctly predicted by the flow-cell model at 5 years was not one of the wells correctly predicted at 50 years but instead was 21% lower than the least-squared DCA regression forecast. The N2 and N3 wells were correctly predicted at 50 years with less than 2% mismatch for both of the wells as compared to the least-squared DCA regression forecasts. However, the 5-year recovery was overpredicted by 17% and 18%, respectively. The N7 well had the largest difference between the two methods with 54% at 5 years and 62% at 50 years. The completion design of the N7 was the closest to the type well with cluster spacing only 20% closer yet the model predicted the biggest difference in production.
The full production rate-time and cumulative production-time plots for each of the wells forecasted with the flow-cell model are shown in Appendix B. As seen with the tabular comparison of Table 12, the flow-cell model tends to overpredict the initial production significantly. The 1-year cumulative production is overpredicted in all of the wells, and is between 20% and 60% higher than the actual production. The consistent, large overprediction of production within the first year will have a large impact on the value metrics (PV10, IRR, PVI, etc.).
What next needs to be answered is this: what are the possible explanations for the discrepancy between the numerically verified flow-cell model predictions and the actual field performance? For example, the flow-cell model relies entirely on the assumption that all of the hydraulic fractures are contributing to the initial rate of production. It is possible that in the actual new wells realized in the field, not all of the hydraulic fractures are equally contributing. This topic will be addressed in more detail in the section below.

Matching Flow-Cell Model Results with the 11 Case Study Wells
With the proof that the numerical reservoir model and the flow-cell model match well over time, the possible causes of the early time overprediction between the flow-cell model forecasts and the actual well performance in the case study wells is analyzed below in more detail.
One of the key assumptions within both the reservoir model and the flow-cell model is that there is a constant bottom hole pressure (BHP). The numerical simulation assumed a bottom hole pressure of 600 psi for the life of the well. The flow-cell model assumes the pressure profile of the well or wells within the type well used to determine the hyperbolic forecast. Gray's correlation, which is used for predominately gas well, was used with the surface casing pressure, gas flow rate, and water flow rate to estimate the bottom hole pressure of all of the case study wells within IHS Harmony software.
A plot of BHP over time is shown in Figure 16. The plot shows a large spread in pressure decline over time. The initial BHP on all of the wells is similar but the pressure decline is significantly different. The pressure decline is similar amongst wells of the same completion design (same colors on the plot) except for one well (N1). The most recent completions (red and purple on the plot) had larger surface facilities and produced without any tubing unlike the rest of the wells. Flowing the wells without any tubing helped lower the BHP much quicker as there was significantly less friction within the production system. The most recent wells (2018/2019) produced below the assumed BHP of 600 psi within 30 days of first production whereas the wells within the type well ranged from 200 days to 1400 days of first production before reaching a BHP of 600 psi. This shows that all of the wells do not have the same BHP over time which may conflict with that assumption of the flow-cell model.

Back Calculation of Flow-Cell Parameters Using Field Production Forecast
The numerical model and flow-cell model both assumed uniform, planar hydraulic fractures that all contribute equally to the production rate. It is probably too idealistic to assume that each hydraulic fracture receives equal water and proppant creating the same fracture half-length and width. As multiple studies have shown, not all fractures are created equal [10,17,22].
To test which proportion of the fractures in the actual wells contributed to production, the inputs of the flow-cell model were treated as unknowns, while the output (actual production forecast) was considered as known. The type well inputs to be matched remained constant, while the number of fractures, fracture spacing, and the fracture halflength in the flow-cell model were varied to calculate the Arps decline parameters (qi, Di, and b-factor) that optimally matched the DCA of the individual wells as shown in the example for one well in Figure 17. Table 13 shows the original inputs to the flow-cell model (as determined from the executed completion design and analytical model match) and the back calculated inputs that match the individual well DCA. Table 14 shows a comparison of the original inputs and the optimal inputs. In all of the cases, the number of fractures were reduced and the fracture spacing was increased when back calculated compared to the original inputs. On average, the number of fractures decreased by 27% and the fracture spacing increased by 38%. The number of fractures is the only parameter used to calculate the initial rate using the flow-cell model equations. In all of the cases with the original number of fractures, which were all assumed to be contributing equally to the actual production, the flow-cell model overestimated the initial rate. This is strong evidence that not all of the fractures planned in the fracture treatment design were created equal, and in fact not all fractures are contributing to the well, according to the proportions given in Table 13. The fraction of non-contributing fracture clusters (non-contributing fractures) ranges between 22% and 36%. This is an important result, which suggests that with the strong emphasis on rapid completions, the well performance is negatively impacted by approximately a quarter of all perforation clusters failing to create fracture systems that would otherwise contribute to production. The conclusion is that past fracture treatment quality assurance is prudent and deserves more attention in There are not any obvious trends between the BHP over time and the wells that had overpredicted or underpredicted production rate or cumulative production. However, if the modern wells had a BHP profile similar to the wells within the type well they would have significantly lower gas rates and lower cumulative production causing the flow-cell model to overpredict early production rate and cumulative production even more than already overpredicted.

Back Calculation of Flow-Cell Parameters Using Field Production Forecast
The numerical model and flow-cell model both assumed uniform, planar hydraulic fractures that all contribute equally to the production rate. It is probably too idealistic to assume that each hydraulic fracture receives equal water and proppant creating the same fracture half-length and width. As multiple studies have shown, not all fractures are created equal [10,17,22].
To test which proportion of the fractures in the actual wells contributed to production, the inputs of the flow-cell model were treated as unknowns, while the output (actual production forecast) was considered as known. The type well inputs to be matched remained constant, while the number of fractures, fracture spacing, and the fracture halflength in the flow-cell model were varied to calculate the Arps decline parameters (q i , D i , and b-factor) that optimally matched the DCA of the individual wells as shown in the example for one well in Figure 17. Table 13 shows the original inputs to the flow-cell model (as determined from the executed completion design and analytical model match) and the back calculated inputs that match the individual well DCA.
In addition, the fracture half-length changed substantially on almost all of the wells from the original inputs to the back calculated inputs. The fracture half-length increased on average 25%. With the amount of fractures reduced, the fracture half-lengths were mostly increased to ensure that the Di and b-factor matched the DCA. The back calculated half-lengths had a much tighter range (548 to 652 ft) than the original input half-length (384 to 776 ft). This suggests that there is too much variability in the history match analytical model in estimating the fracture half-length for the flow-cell model. Figure 17. Regression of the monthly gas rate and the cumulative gas volumes to back calculate the inputs of the flow-cell model for an optimal match of N4 (without any minimum decline).   Figure 17. Regression of the monthly gas rate and the cumulative gas volumes to back calculate the inputs of the flow-cell model for an optimal match of N4 (without any minimum decline).  N7  61  776  81  48  560  104  N4  133  384  58  96  548  80  N5  159  652  47  122  640  61  N6  162  439  46  123  641  61  N2  176  445  36  113  652  56  N3  161  446  40  109  645  59   Table 14 shows a comparison of the original inputs and the optimal inputs. In all of the cases, the number of fractures were reduced and the fracture spacing was increased when back calculated compared to the original inputs. On average, the number of fractures decreased by 27% and the fracture spacing increased by 38%. The number of fractures is the only parameter used to calculate the initial rate using the flow-cell model equations. In all of the cases with the original number of fractures, which were all assumed to be contributing equally to the actual production, the flow-cell model overestimated the initial rate. This is strong evidence that not all of the fractures planned in the fracture treatment design were created equal, and in fact not all fractures are contributing to the well, according to the proportions given in Table 13. The fraction of non-contributing fracture clusters (non-contributing fractures) ranges between 22% and 36%. This is an important result, which suggests that with the strong emphasis on rapid completions, the well performance is negatively impacted by approximately a quarter of all perforation clusters failing to create fracture systems that would otherwise contribute to production. The conclusion is that past fracture treatment quality assurance is prudent and deserves more attention in future field completions.
In addition, the fracture half-length changed substantially on almost all of the wells from the original inputs to the back calculated inputs. The fracture half-length increased on average 25%. With the amount of fractures reduced, the fracture half-lengths were mostly increased to ensure that the D i and b-factor matched the DCA. The back calculated half-lengths had a much tighter range (548 to 652 ft) than the original input half-length (384 to 776 ft). This suggests that there is too much variability in the history match analytical model in estimating the fracture half-length for the flow-cell model.

Additional Explanations
Both the numerical model and flow-cell model assume uniform, planar hydraulic fractures that do not connect to any other hydraulic fractures or any natural fractures. While there is not any evidence of natural fractures within this particular lease region, it is plausible that there is interwell communication through connected hydraulic fractures. Within the dataset, there is extra produced water on existing wells from offsetting completions meaning the stimulated rock volumes are connected.
For example, there was extra produced water on existing wells from offsetting completions. For example, the N1 had a peak of water production of 1000 bbls/day at initial production but had decreased to an average of 15 bbls/day after 4 years of production. After the offset N2 completion, the water production on the N1 spiked to 380 bbls/day and continued to be elevated for the next 3 months. Unfortunately, neither well had a shut-in test to monitor pressures and rates during interwell communication. There is not clear evidence of connected fractures after the initial completion meaning it is unknown if the drained rock volumes were affected.
The fracture half-length for the best match of the numerical simulation was 600 ft (Table 8) while the interwell spacing is 1300 ft or a half-length of 650 ft. The original inputs used in the flow-cell model had an average fracture half-length of 520 ft (Table 13), while the back calculated inputs have a significantly higher average fracture half-length of 615 ft (Table 13). Clearly, to assume that each hydraulic fracture has the same fracture half-length may be too idealistic.
For example, fractures from neighboring wells may be interconnected, particularly since the sum of two fracture half-lengths almost equals the 1300 ft well spacing. Chin et al. [23] showed in his study that the reservoir model for recent wells with modern completion designs did not work until connected fracture networks were incorporated into the reservoir model. The interwell communication through connected hydraulic fractures should affect the decline rate and the b-factor by increasing the initial decline and decreasing the b-factor compared to a well without connected hydraulic fractures, which will decrease overall recovery.

Verifying the Flow-Cell Model with a Reservoir Model
A CMG numerical simulation was built that successfully modeled the production of the chosen type well. The history match allowed a few different reservoir properties and fracture properties to vary to accurately model the production. Most of the properties did not change much between the initial guess and the best match. For example, permeability only changed from 100 to 130 nD while matrix porosity only change from 7 to 6.98%. Both were allowed to vary significantly ( Table 8). The small change means that the measurements from the nearby core sample were accurate. The fracture width (0.000833 ft to 0.0002255 ft) and fracture permeability (20,000 mD to 5655 mD) changed significantly from the initial guess to the best match ( Table 8). The initial guess was based on the assumption that the fractures were infinitely conductive meaning that the pressure drop along the fracture would be approximately zero (F cd > 30). However, both the fracture width and fracture permeability decreased by a factor of 4 causing the dimensionless fracture conductivity to drop to 16.6, as shown in Table 15. With a dimensionless fracture conductivity of 16.6, the fractures in the numerical model are no longer infinitely conductive. One of the main assumptions of the flow-cell model is that there is near infinite fracture conductivity because fluid flow into the fracture is assumed to be equal along the entire length of the fracture. Using the modified Pratt Equation (9), a dimensionless fracture conductivity of 30 (infinite conductivity) gives an effective fracture half-length that is 95% of the flowing fracture half-length [26]. A dimensionless fracture conductivity of 16.6 gives an effective fracture half-length that is 91.5% of the flowing fracture half-length, which is close to infinite conductivity. The flow-cell model assumption of near infinite fracture conductivity is therefore still met. x Separately, non-infinitely conductive fractures may occur in the study wells, due to a less than optimal completion procedure. Fracture conductivity controls early time transient production. As all of the wells are still within transient flow (specifically linear flow as seen in Section 3.4) after 4 years of production, it becomes clear that transient production is important to the production rate. Increasing the fracture conductivity should be investigated. Using too small of sand to prop the fractures, poor sand quality, too little water pumped, too much friction when pumping, too large of a leak-off-rate are some of the potential causes of fracture treatment procedure delivering a fracture conductivity that is lower than intended.

Verifying the Reservoir Model: Optimal Match of Numerical Simulation
The original comparison of the numerical simulation and the type well was given in Figure 12. The original comparison of the data was plotted on a logarithmic scale which is great for viewing data across different orders of magnitude. However, logarithmic plots have a tendency to suppress the mismatch, because of the vertical scale compression. The numerical simulation compared to the type well forecast had an average absolute relative error of 2.3% and 15.2% over 30 years for cumulative gas production and daily gas production, respectively. The percentage of mismatch over time of the numerical simulation against the type well forecast is shown in Figure 18 and Table 16. Figure 12. The original comparison of the data was plotted on a logarithmic scale which is great for viewing data across different orders of magnitude. However, logarithmic plots have a tendency to suppress the mismatch, because of the vertical scale compression. The numerical simulation compared to the type well forecast had an average absolute relative error of 2.3% and 15.2% over 30 years for cumulative gas production and daily gas production, respectively. The percentage of mismatch over time of the numerical simulation against the type well forecast is shown in Figure 18 and Table 16. Figure 18. The average error between the numerical simulation and the type well for gas production rate and cumulative production. Figure 18. The average error between the numerical simulation and the type well for gas production rate and cumulative production. Initially, the numerical simulation had a lower forecasted gas rate than the type well, which lasted for the first year. The numerical simulation forecast remained higher until 14 years of production. The gas rate of the numerical simulation finished lower at 30 years of production. The largest mismatch of the gas rate was at 30 years of production where there was a 23% difference between the numerical simulation and the type well forecast. The average absolute mismatch for the first five years of production was the same as the 30-year average of 14.8%.
Due to the lower initial gas rate, the numerical simulation had a lower cumulative gas production for the first 6 years than the type well forecast. The numerical simulator then remained with a higher cumulative gas production until 22 years of production when the type well forecast surpassed it. The numerical simulation cumulative gas production at 30 years ended lower than the type well forecast (by 2.3%). The largest mismatch on the cumulative production occurred after 0.75 years of production where there was a 13% difference. The average mismatch for the first five years of production was 6.6%.
While the gas rate mismatch is significant, the cumulative production curve was extremely accurate. While an individual month's production may be off, the total production and therefore the total revenue and total cash flow will be correct. There should be a little caution with any present worth calculations done with this forecast; however, the total estimated recovery is accurate.

Flow-Cell Model
The flow-cell model was tested against the numerical reservoir simulator (in Section 4.3) for use in dry gas wells with four different fracture spacing. Although the overall matches were satisfactory (Figure 15), the flow-cell model has a tendency to overpredict the initial production compared to the numerical simulation especially within the first 5 years (Figure 15). On the 150 ft fracture spacing test, the flow-cell model did not overpredict early production. For the 50, 70, and 100 ft fracture spacing models, the flow-cell model overpredicted early production. However, it should be noted that the numerical simulator tends to underpredict the early production ( Figure 18).
There appears to exist a correlation between the flow-cell overprediction and fracture spacing: the lower the fracture spacing the higher the overprediction in the early time production data; the tighter the fracture spacing the more drastic and sooner the overprediction occurred (Figure 15). This early overprediction of the well rates, when combined with an economic appraisal, may possibly cause slightly inflated rates of returns due slower discounting in the present value of money.
Another source of any observed mismatch between flow cell-model and real well performance could be due to an earlier mismatch between the history matched decline curve for the type well and the actual type well (due to inclusion of spurious early data points that may need to be omitted to get DCA parameters giving a better long term fit).
Coupling of the flow cell model with an economic module can help find the well spacing and fracture treatment parameters that lead to the economic optimum (NPV).

Conclusions
This study started off by validating the flow-cell model with a numerical simulation of a Marcellus gas well. The numerical simulator was successfully matched to the historical performance of one of the wells. The fracture spacing of both the flow-cell model and the numerical simulator were varied from 50 to 150 ft. The flow-cell model and numerical simulator closely matched on the gas flow rate and cumulative gas production over time for all of the various fracture spacing. Compared to the numerical simulation, the flow-cell model proved to be a quick and accurate model.
The flow-cell model was then used in a case study of 11 horizontal, multi-fractured, dry gas Marcellus wells. Of the 11 wells, 5 wells with similar completion designs were utilized to create a type well. The flow-cell model used the type well to predict the gas flow rate and the cumulative production for the remaining 6 wells with different completion designs. The flow-cell model predictions were compared to the historical production of the 6 wells.
Overall, the flow-cell model did not match the gas flow rate or cumulative production of actual well results in the case study. The flow-cell model tended to overpredict the early time production benefits of increasing the number of fractures. As validated by the numerical reservoir model, the flow-cell model is an accurate reservoir model. The mismatch of the flow-cell model predictions and the performance of the case study wells must come from some of the model assumptions or some of the inaccuracies in the inputs. A few sources of mismatch between the two data sets were identified.
The study initially assumed that 100% of the perforation clusters created fractures in the field wells. However, in the literature, other operators in the Appalachian basin assumed between 33% and 60% of the perforation cluster created fractures. The flow-cell model was used in reverse, where the model was given the forecast for the new well to see the fracture dimensions needed to create that forecast. In all of the wells, the number of fractures decreased (on average~27%), and the fracture spacing increased (on average 38%). As the number of fractures is the only parameter to calculate the initial rate using the flow-cell model, the early time overprediction is explained by having less than 100% cluster efficiency in the field wells.
Within the flow-cell model, the fracture half-length is also important for predicting the size of the box drained to achieve a particular EUR. In the use of the flow-cell model to determine the fracture dimensions based on the forecast, all of the fracture half-lengths changed and became closer to the type well averages for the fracture half-length. The IHS harmony model used to estimate the fracture half-length does not have enough uniqueness to provide accurate fracture half-lengths, which original inputs may have caused imprecise forecasts and EURs from the flow-cell model.
Lastly, one of the key assumptions within the flow-cell model is that there is a constant bottom hole pressure (BHP). The initial BHPs on all of the 11 case study wells were similar, but their pressure declined over time, and in different ways due to different production systems being used. The variations in the pressure decline over time did not exhibit any obvious trends between which wells had overpredicted or underpredicted performance; however, the variable pressure decline may be another source of mismatch between the flow-cell model results and the actual well performance.
Finally, this study lists the following insights: • The flow-cell model has been validated for a horizontal, multi-fracture, gas well based on the close match to the numerical simulation well.

•
The flow-cell model forecasted rates mismatched performance of actual wells in the case study.

•
The inaccuracies of the flow-cell model forecasts, practically the early time overprediction in the case study, must come from some of the assumptions or some of the inputs. The flow-cell model matches the actual well performance when about a quarter of the perforation clusters were assumed inactive.

•
In addition, estimates of fracture half-lengths used in the flow-cell model were adjusted to achieve better matches with the actual well performance. The adjusted half-lengths had a significantly tighter range (100 ft) than those estimated by the IHS Harmony model (400 ft).

•
The large variance in pressure decline over time between the wells may be another source of the mismatch between flow-cell model forecasts and actual well performance.

•
The mismatch earlier observed by the operator between the numerical reservoir model forecast and actual well production may have been due to the failure of perforation clusters to perform.

•
In this respect, the flow-cell model provided flexible/versatile enough for reverse use and history match of the actual well data to reveal which completion parameters and in what proportion needed to be changed to achieve good matches. EURs. An engineering judgment was therefore made to limit the b-factor between 1.0 and 1.13. An example of the regression fit on the monthly production data for N1 is shown in Figure A1. The Arps decline parameters for the fits for all of the wells studied are shown in Table A1. The more recent wells have significantly higher flow rates, normalized flow rates, and decline rates. Notably, 10 out of the 11 wells share the same b-factor of 1.13.    Table A2. The 5-year recovery is important, since most of the value of the well comes from production within this time frame while the 50-year recovery represents the reserves estimates associated with each well. A cumulative distribution frequency of the normalized recovery is shown in Figure A2. All of five of the modern completion designs (completed after 2015) were in the 50th percentile or better for the 5-year recovery volume; however, only three out of the five wells were in the 50th percentile or better for the 50-year recovery. This means that the new completion techniques are accelerating production but are not necessary increasing the EUR.  Figure A2. Cumulative distribution frequency vs. normalized recovery over 5 years (a) and 50 years (b).

A.3. Type Well Construction
An empirical type well was generated from the five wells with the same 2015 completion design (5 clusters at 100 ft spacing in one stage with 2000 lbs of proppant per foot of lateral) to use in the flow-cell model for estimating the production of all of the others wells. The production rate of the empirical type well was plotted with the individual wells ( Figure A3). The decline parameters for the type well are shown in Table A3. To ensure that the type well is indeed representative of the early time and the late-time behavior of the wells, the type well should be close to the average, normalized 5-year recovery and 50-year recovery of the individual wells. Table A4 shows that this is the case, with less than 1% error in both time frames meaning the type well is representative of the average well with 5 clusters at 100 ft spacing in one stage with 2000 lbs of proppant per lateral foot.

Appendix A.3. Type Well Construction
An empirical type well was generated from the five wells with the same 2015 completion design (5 clusters at 100 ft spacing in one stage with 2000 lbs of proppant per foot of lateral) to use in the flow-cell model for estimating the production of all of the others wells. The production rate of the empirical type well was plotted with the individual wells ( Figure A3). The decline parameters for the type well are shown in Table A3. To ensure that the type well is indeed representative of the early time and the late-time behavior of the wells, the type well should be close to the average, normalized 5-year recovery and 50-year recovery of the individual wells. Table A4 shows that this is the case, with less than 1% error in both time frames meaning the type well is representative of the average well with 5 clusters at 100 ft spacing in one stage with 2000 lbs of proppant per lateral foot.
( Figure A3). The decline parameters for the type well are shown in Table A3. To ensure that the type well is indeed representative of the early time and the late-time behavior of the wells, the type well should be close to the average, normalized 5-year recovery and 50-year recovery of the individual wells. Table A4 shows that this is the case, with less than 1% error in both time frames meaning the type well is representative of the average well with 5 clusters at 100 ft spacing in one stage with 2000 lbs of proppant per lateral foot.    The actual production, least-squared regression forecast, and the flow-cell model forecasts for all six wells with different completion designs than the type well are shown in Figures A4-A9 (using monthly production vs. time and cumulative production vs. time plots).

Appendix C. Original Gas in Place and Recovery Factors
The estimated ultimate recovery (EUR) of the case study wells was not tied to the recovery factors when production forecasting with the decline curve analysis (Appendices A and B). It is always good to compare the EURs to the original gas in place (OGIP) and determine the recovery factors to avoid unrealistic EURs.
The OGIP is a volumetric calculation of how much gas exists within the reservoir. There are two parts to the calculation: free and adsorbed. Free gas occupies the pore space whereas the adsorbed gas is a thin film of gas attached to the surface area of the matrix. The adsorbed gas needs a low pressure to be able to desorb from the surface of the matrix and does not typically contribute to the production in the early times. The equations in field units for the free, adsorbed, and total OGIP are shown in Equations (A1), (A2) and (A3), respectively.

Appendix C. Original Gas in Place and Recovery Factors
The estimated ultimate recovery (EUR) of the case study wells was not tied to the recovery factors when production forecasting with the decline curve analysis (Appendices A and B). It is always good to compare the EURs to the original gas in place (OGIP) and determine the recovery factors to avoid unrealistic EURs.
The OGIP is a volumetric calculation of how much gas exists within the reservoir. There are two parts to the calculation: free and adsorbed. Free gas occupies the pore space whereas the adsorbed gas is a thin film of gas attached to the surface area of the matrix. The adsorbed gas needs a low pressure to be able to desorb from the surface of the matrix and does not typically contribute to the production in the early times. The equations in OGIP _Total = OGIP _Free + OGIP _Adsorbed (A3) Within those equations, A is area, H is pay zone height, S gi is the initial gas saturation, φ is the porosity, B gi is the initial gas formation volume factor, ρ b is the bulk density, V l is the Langmuir volume, p i is the initial pressure, and p l is the Langmuir pressure.
The total OGIP was calculated for each of the wells based on the individual completed lateral length and the interwell spacing. The assumed Langmuir pressure and Langmuir volume come from a paper about adsorption from the Marcellus from the National Energy Technology Laboratory [28]. All of the other inputs for the free and adsorbed OGIP calculations can been seen in Table A7. The free gas OGIP makes up~40% of the total OGIP. The total OGIP compared to the 50-year EUR from the Arps DCA and the flow model are shown in Table A8. The key assumption is that the well will drain all of the volume within the evenly split lateral well spacing. The range of the recovery factors for the Arps DCA is 27% to 46% with an average recovery factor of 38%. A conventional gas reservoir with gas expansion as the reservoir drive typically recover 80-90% of the OGIP [29]. The Due to significantly lower permeability among many other differences, unconventional reservoirs are not expected to achieve similar recovery factors as conventional reservoirs; therefore, these recovery factors are reasonable.
The range of recovery factors for the flow-cell model is 29% to 64% with an average recovery factor of 42%. The flow-cell model assumes when the fracture half-length increases in the new well from the type well the EUR increases too but the recovery factor stays constant. The two wells, N7 and N5, that had fracture half-lengths higher than the type well also had the largest recovery factors with the flow-cell model. The input fracture half-length of the flow-cell model is important for determining the EUR of the new well.
The N7 and the N5 had an inferred fracture half-length of 776 ft and 652 ft, respectively (see Table A6, Appendix A), compared to the type well average fracture half-length of 509 ft. The lateral well spacing for those two wells is 1300 ft, meaning that the inferred fracture half-lengths are extending beyond the 1300 ft lateral spacing (also used in the OGIP calculation). This disparity suggests that the fracture half-lengths are too high. The other four wells with lower fracture half-lengths than the type well average had lower recovery factors with the flow-cell model between 29% and 33% compared to the type well average of 37%.