Forecast of economic gas production in the Marcellus Shale

Abstract

We apply a hybrid data-driven and physics-based method to predict the most likely futures of gas production from the largest mudrock formation in North America, the Marcellus Shale play. We first divide the ≥100,000 mi2 of the Marcellus Shale into four regions with different reservoir qualities: the northeastern and southwestern cores and the noncore and outer areas. Second, we define four temporal well cohorts per region, with the well completion dates that reflect modern completion methods. Third, for each cohort, we use generalized extreme value statistics to obtain historical well prototypes of average gas production. Fourth, cumulative production from each well prototype is matched with a physics-based scaling model and extrapolated for two more decades. The resulting well prototypes are exceptionally robust. If we replace production rates from all of the wells in a given cohort with their corresponding well prototype, time shift the prototype well according to the date of first production from each well, and sum up the production, then this summation matches rather remarkably the historical gas field rate. The summation of production from the existing wells yields a base or do-nothing forecast. Fifth, we schedule the likely future drilling programs to forecast infill scenarios. The Marcellus Shale is predicted to produce 85 trillion SCF (TSCF) of gas from 12,406 existing wells. By drilling ∼3700 and ∼7800 new wells in the core and noncore areas, the estimated ultimate recovery is poised to increase to ∼180 TSCF. In contrast to data from the Energy Information Administration, we show that drilling in the Marcellus outer area is uneconomic.Natural gas plays an essential role in the possible transitions to clean energy. Today, natural gas fulfills one-fourth of the global primary power demand (BP, 2020), and its importance to the global power supply mix is predicted to only increase in the next two decades. In the United States, natural gas provides one-third of the total primary power demand (Energy Information Administration, 2020b). In mid-2020, the United States produced nearly 110 BCF of natural gas per day, which is almost twice the production rate of 15 yr ago (Enverus, 2021). This significant increase in natural gas output in the United States has only been possible with the “unconventional resource revolution” over the last two decades, during which operators have learned how to produce gas from the extremely low-permeability—unconventional—mudrock or so-called shale formations by advancing horizontal drilling and hydraulic fracturing technologies. Today, the eight major shale plays in the United States are responsible for nearly 70% of the total natural gas output. The Marcellus Shale, the largest gas shale play in North America, contributes one-third of the total United States shale gas production, producing more than 25 BCF of natural gas per day. As of the writing of this paper, the Marcellus has produced 50 trillion SCF (TSCF) of natural gas, equivalent to 8.3 billion bbl of oil.Despite the key role of the Marcellus in securing national power demand, a significant disparity exists in the predictions of the total amount of technically recoverable gas contained in the Marcellus. A study by the Bureau of Economic Geology (BEG) at The University of Texas at Austin (Ikonnikova et al., 2018), predicted the technically recoverable resources (TRRs) of the Marcellus to be 560 TSCF. The US Department of Energy’s Energy Information Administration (EIA) recently published a prediction of 310.6 TSCF of technically recoverable gas (Energy Information Administration, 2020a). Lastly, the US Geological Survey (USGS) estimated only 96.5 TSCF of total undiscovered resources in the Marcellus. This divergence of total producible resource predictions in the Marcellus shows just how difficult and uncertain the resource assessment in unconventional reservoirs is. Conventional reservoirs can be evaluated by identifying their structures, traps, and water–gas contacts. Assessment of recoverable resources in unconventional reservoirs is highly uncertain and unfeasible with use of the traditional volumetric methods (Schmoker, 2002; Haskett and Brown, 2005; Cipolla et al., 2011; Zou, 2017). In unconventional reservoirs, hydrocarbon accumulations exist continuously and span huge areas of shale or tight carbonate–sandstone formations. However, because of the extremely limited reservoir connectivity, thousands of horizontal hydraulically fractured wells must be drilled to achieve a desirable rate of production. Therefore, the risk of hydrocarbon finding, which is key to assessing conventional reservoirs, becomes insignificant relative to the risk of commerciality (Fulford, 2016).In many unconventional plays, better-producing zones are concentrated in smaller geographic and geologic zones that are often called “sweet spots” (Hughes, 2013). Another factor that compounds the uncertainty of unconventional resource assessment is the advancement of completion technologies (Weijermars, 2015). Modern wells with new completions generally outperform vintage wells in the same reservoir, but with traditional completion systems. These new wells are completed with intensive hydraulic fracturing treatments (more of larger hydraulic fractures).In this paper, we pursue our tested integrated approach, which minimizes the uncertainty of fieldwide resource assessment in the Marcellus. Following USGS (Higley et al., 2019), we focus only on the Middle Devonian Marcellus layer. There may be additional potential production coming from adjacent formations, such as the Upper Devonian Shale. However, as of January 2021, there were only 399 wells producing from the Upper Devonian that produced negligible gas in comparison with the 15,000+ wells drilled in the Middle Devonian Marcellus (see, for example, Milici and Swezey, 2014, for a full stratigraphic cross section/type section of the Appalachian Basin). By including fundamentals of shale geology, advancement of well completion technologies, physics of natural gas production from the hydraulically fractured shales, and economics of drilling projects, our approach is different from the volumetric approach of USGS (Higley et al., 2019). According to the classification of the resource assessment approaches for unconventional formations (Schmoker, 2002), our method belongs to his second category that relies heavily on production performance of existing wells, instead of the highly uncertain volumetric estimates in category one.First, we divide the Marcellus Shale play into several spatiotemporal well cohorts in which gas production is statistically uniform. Second, for each cohort, we construct a hybrid data-driven and physics-based well prototype based on the generalized extreme value (GEV) statistics (Fréchet, 1927; Weibull et al., 1951; Gumbel, 1954) and on physical scaling (Patzek et al., 2013, 2014, 2017; Eftekhari et al., 2020; Saputra et al., 2020, 2021a; Saputra, 2021). We have previously introduced the hybrid GEV + physical scaling method in Patzek (2019), Patzek et al. (2019), and Saputra et al. (2019), in which we constructed well prototypes for the Barnett Formation and Bakken Formation shales. To check whether we have the robust well prototypes that best represent at times more than 1000 wells per cohort, we replace production from each existing well with its corresponding well prototype. Since the summation of production from the time-shifted well prototypes matches the total historical field rate, the extrapolated future production serves as our base or do-nothing forecast, which gives a robust estimated ultimate recovery (EUR) from all of the existing wells in the Marcellus.To obtain the fieldwide EUR, we also calculate the most plausible numbers of wells that can be drilled in each subregion of the Marcellus, and schedule future drilling programs. Replacing each new well drilled with a prototype appropriate for a subregion, we finally obtain the robust infill forecasts. Finally, we perform a full economic analysis of each infill scenario, and determine which subregions of the Marcellus are technically and economically acceptable and which are not. Our predicted TRR value, 178 to 200 TSCF, is similar to that of USGS (Higley et al., 2019), yet it is much lower than the predictions by EIA (Energy Information Administration, 2020a), 310.6 TSCF, and BEG (Ikonnikova et al., 2018), 560 TSCF. Here, we show how we reached this conclusion.In 1839, James Hall, from the New York State Geological Society, discovered a black shale outcrop in the town of Marcellus, Onondaga County, New York, and named it the Marcellus Shale. In the 1930s, some vertical wells penetrated the black shales of the Marcellus Shale at approximately 1000 ft deep and encountered a strong influx of gases. However, these flows came only from isolated “gas pockets” and could not be sustained (Harper, 2008; Carter et al., 2011). In the late 1970s and early 1980s, more exploratory wells were drilled and cored to map the distribution of the Marcellus in Pennsylvania. The Pennsylvania Geological Survey concluded that the Marcellus has excellent potential as an important gas reservoir despite its extremely low permeability (Harper, 2008). Up to 2005, only 100 vertical wells had been drilled and produced less than 1.5 million SCF (MSCF)/day. A few years after, led by Range Resources, the new horizontal hydraulically fractured wells that made the Barnett Shale so successful were applied massively in the Marcellus. Today, there are 15,000+ horizontal Marcellus wells producing 25 billion SCF (BSCF)/day, nearly one-third of United States natural gas output.In addition to being the largest gas producer in the United States, the Marcellus is the largest United States shale gas play by area. Figure 1A, published by EIA (Popova, 2017) shows the area of the Marcellus, approximately 104,000 mi2, which extends over several states—nearly the entirety of West Virginia, three-fourths of Pennsylvania, and approximately half of New York and Ohio. This black Marcellus Shale is among the richest mudrock plays in the United States, with total organic carbon of up to 20 wt. %. It was deposited in the Appalachian foreland basin during the mid-Devonian, circa 385 Ma (see Table 1 for a simplified stratigraphic column of the Appalachian Basin [Milici and Swezey, 2014]).The EIA outlines 58,000 mi2, depicted with the thick dashed line in Figure 1A to delineate the Marcellus play extent, where drilling should be profitable to the operators. To show depth and thickness variations in the Marcellus Shale, we have constructed three cross sections—AA′, BB′, and CC′—based on the structure and isopach maps of the Marcellus from Popova (2017) (see Figure 1B). Based on the absolute thickness and the absolute depth of the Marcellus, we can roughly locate two sweet spots, one in northeastern Pennsylvania and the other at the border between northwestern Virginia and southwestern Pennsylvania. The maps of vitrinite reflectance (Ro) and hydrogen index (HI) from Dolson (2016) also confirm that these two sweet spots are thermally mature for producing dry gas and natural gas liquids, with the Ro values of 1.5% to 2.5% and HI values of 4 to 52.To assess technically recoverable resources in the Marcellus, USGS (Higley et al., 2019) defined six geology-based assessment units (AUs) shown in Figure 2. According to USGS, the potential extent of the Marcellus Shale resources is based on the minimum depth of 1000 ft, minimum thickness of 25 ft, and minimum Ro of 0.5%. The total area of the six AUs is almost identical to the area of the Marcellus Shale by EIA in Figure 1A. The best AUs are denoted by the so-called interior affixes, and their total area is much smaller than that of the EIA’s play extent in Figure 1A. Also note that the two sweet spots identified in Figure 1 are, respectively, inside of the northern interior AU and the southern + southwestern interior AUs.In this paper, we set out to construct several AUs similar to those of USGS (Higley et al., 2019). We base our units solely on well quality, captured as well cohorts. To account for the variations of reservoir quality (geology) in the Marcellus, but also for the changes in completion technologies over time, we make these well cohorts spatiotemporal. For every well in each cohort, the goal is to have similar gas production histories so that each well sample is statistically uniform. The main data sources for this study are obtained from Enverus (2021) and FracFocus (2021), and detailed in the Appendix, Table 9.First, we divide the Marcellus play extent in space into four reservoir quality areas: the northeast core, southwest core, Marcellus noncore, and Marcellus outer. The two core areas are sweet spots, defined as the envelopes of wells with normalized cumulative production of at least 80 MSCF (see Figure 3). To make a fair comparison between the old and young wells, we divide their cumulative production by the number of elapsed months on production. The shape of the northeast core area follows nicely the USGS northern interior AU in Figure 2, whereas the southwest core area is located in the southern + southwestern interior AUs. These comparisons demonstrate that our core areas also capture the essential geology of the Marcellus, as in Higley et al. (2019). The Marcellus noncore area encompasses the mid-tier wells and is defined as the difference between the effective area and the two core areas. The Marcellus outer area is further defined as the difference between the EIA play boundary and the effective area. It contains subpar wells with normalized cumulative production of less than 20 MSCF.Second, we subdivide the well cohorts in the Marcellus into four completion date intervals: 2009–2012, 2013–2014, 2015–2016, and 2017–2020. This temporal division is based on Figure 4, which shows 10 yr of advancements in well completion technologies in the Marcellus Shale play. From Enverus, we obtain the lateral length data (feet) for 12,401 wells, plotted as gray dots in Figure 4A (Enverus, 2021). Using GEV statistics, we obtain the annual P10 (upper bound), mean (expected value), and P90 (lower bound). From FracFocus, we obtain only 9926 data points for fracturing water volume (gallons) and mass of proppant (pounds) (FracFocus, 2021). We further divide the water volumes and proppant masses for each well in FracFocus by their corresponding lateral lengths from Enverus and obtain the fracturing water intensity (gallons per foot) and proppant intensity (pounds per foot) in the Marcellus (see Figure 4B, C). Since the number of fracture stages for each well is not available, we approximate it by dividing the lateral lengths from Figure 4A with the hydraulic fracture spacing reported from the literature (Gerdom et al., 2013; Eclipse Resources, 2017; Walzel, 2019). In Figure 4D, along with the calculated number of fracture stages, we also show the true values from the literature. From the four time series in Figure 4, we observe a massive completion technology advancement in the Marcellus. Over a decade, the operators almost doubled the lateral lengths and quintupled the number of fracture stages. The increases in fracturing water intensity and proppant intensity may also indicate that more fractures per unit area have been created; hence, reservoir exposure and drainage radius per well have increased.In summary, we divide the Marcellus Shale play into 16 spatiotemporal well cohorts in total (see Figure 5). Table 2 details the sample size numbers included in the cohorts. It is obvious that over the last decade, operators have drilled only in the sweet spots in the two small core areas. Only 12% of the total horizontal wells have been drilled in the Marcellus noncore, and <1% have been drilled in the Marcellus outer.After dividing the Marcellus into 16 spatiotemporal well cohorts, we construct historical well prototypes for each cohort using the GEV statistics. The GEV distributions are used to model numerous extreme events in nature, such as extreme rainfall over a long period, floods, earthquakes, appearance of large pores or fractures in a rock mass, or production of hydrocarbons from fractured shales (Patzek et al., 2019). In statistics, according to the extremal types theorem, the only stable distribution to model such extreme events is GEV distribution, which is able to capture the maxima of long (finite) sequences of random variables. Because measurable gas production in a shale well requires volumetric fracturing of reservoir rock or complex fractures, each production record is an extreme event (Patzek, 2019; Patzek et al., 2019). Therefore, the distributions of annual gas production from the hydrofractured horizontal wells in the Marcellus are accurately fit with GEV distributions.The previous research by Patzek et al. (2019) and Saputra et al. (2019) shows the superiority of GEV distribution in matching annual gas and oil production in the Barnett and Bakken Shales, compared with the widely used lognormal distribution. The GEV distribution is a generalization of the three extreme value distributions (Fréchet, 1927; Weibull et al., 1951; Gumbel, 1954), each depending on three parameters. In contrast, a normal or lognormal distribution only has two parameters, μ and σ⁠. In a symmetrical distribution, μ is equal to the mean and σ is the standard deviation. A nonsymmetrical distribution such as GEV may also depend on a third parameter. In a GEV distribution the location parameter is μ⁠, the scale parameter is σ⁠, and the shape parameter is ξ⁠. The value of ξ determines which of the three extreme distributions results—Fréchet (⁠ξ>0⁠), Weibull (⁠ξ<0⁠), or Gumbel (⁠ξ=0⁠). In most shale production cases, Fréchet is a better fit because we find that the shape factor, ξ⁠, of annual oil or gas production distributions in several major shale plays in the United States is positive and between 0.1 and 0.2.Note that a good fit of data with a GEV distribution usually has ξ << 1/2, because variance becomes infinite for ξ ≥ 1/2 (see Patzek et al., 2019; Saputra et al., 2019).Figure 6 demonstrates how we obtain a historical well prototype for each of the 16 Marcellus well cohorts in Table 2. For simplicity, this example focuses on cohort-7, which is the southwest core (2015–2016), with 1147 wells. First, we construct subcohort-7-j, which contains all of the wells in cohort-7 with at least j years on production. For example, there are 1145 wells in the subcohort-7-1 because the other 2 wells in cohort-7 are in the early phases of production, with less than 1 × 12 months on production. Similarly, the subcohort-7-5 has only 450 wells because the other 697 wells in cohort-7 are less than 5 yr old. No subcohort-7-6 is found because cohort-7 only contains wells that were completed in 2015 to 2016.Second, for each subcohort, we plot the histogram of annual gas rates (BSCF/yr). We then fit a GEV distribution to this histogram, using the GEV PDF in equation 1 with three fitting parameters: μ⁠, σ⁠, and ξ⁠. We also display the lognormal distribution fit to show that GEV is better fit to the widely used lognormal distribution in matching annual gas production in shales. The second row in Figure 6 shows the contour plots of μ and σ with a corresponding 95% confidence interval. The fewer wells we have in a sample, the wider the uncertainty contour. Because of the fit uncertainty, we discard the subcohorts with fewer than 40 wells. The third row in Figure 6 shows the CDF that can be calculated using equation 2. From the CDF, we can quickly project the P10 (upper bound, or 90th percentile), P50 (median), and P90 (lower bound, or 10th percentile). Using equation 3, we can calculate the mean of each subcohort. Notice that the calculated mean is always higher than the median if we have a positive skewness in the distribution (the heavy tail is on the right-hand side).Third, we plot the means, P10⁠, and P90 of the annual gas rate (BSCF/yr/well) of each subcohort versus years on production and calculate the cumulative production (see Figure 6). At this point, we have a complete historical well prototype for cohort-7. We repeat these steps for the other 15 well cohorts in the Marcellus. Each well prototype is constructed solely from the historical production data of existing wells, and it is limited to 5 yr for the cohort-7 wells. In the next step, we use a physical scaling approach to extend the well prototypes for up to two or three more decades. (See SOM-1, supplementary material available as AAPG Datashare 173 at www.aapg.org/datashare, for the individual plots of GEV PDF, GEV contour, and GEV CDF for all 16 cohorts of the Marcellus.) Our approach to obtaining well prototypes from the expected values of GEV distributions is highly predictive. In Figure 7, we demonstrate a cross-validation that GEV means are remarkably stable, even if we omit 50% of the data.In addition to the robust GEV statistics, we incorporate physical scaling to extend the historical GEV mean well prototypes for two or three more decades (Patzek et al., 2013, 2014, 2017, 2019; Eftekhari et al., 2018, 2020; Saputra et al., 2020). We adopt this physics-based method instead of the purely empirical industry standard hyperbolic decline curve analysis method (Arps et al., 1945), which tends to overestimate future production (Olson et al., 2019; Haider et al., 2020; Saputra et al., 2020).The physical-scaling method we use here was originally developed by Patzek et al. (2013, 2014) as a model of nonlinear diffusion of gas pseudo-pressure during gas production from the hydraulically fractured horizontal wells in the Barnett Shale (see Figure 8). They have established that gas production from ≥8000 wells in the Barnett follows a monotonic master curve, along which cumulative gas production grows initially as a straight line versus the square root of time on production (transient period). Later, after 5 to 10 yr, this cumulative production bends down due to an exponential decrease in production rate in pseudo-steady-state flow. This master curve was obtained by solving numerically a nonlinear pseudo-pressure diffusion equation between two consecutive hydraulic fracturing stages. Later, the research by Patzek et al. (2017) and Saputra et al. (2020) produced an exact analytical solution of the physical scaling method for the horizontal black-oil wells in the Eagle Ford and Bakken Shales. We demonstrated that the physical-scaling method is significantly faster than full-scale numerical reservoir simulation, yet it has predictive power that is similar to numerical simulation and other more complex analytical solutions (Saputra and Albinali, 2018).To account for late-time pressure diffusion and gas production, exterior flow from outside the stimulated reservoir volume (SRV) was added to the master curve. This flow appears as an additional square-root-of-time flow scaled by the ratio of hydrofracture spacing to fracture length (Eftekhari et al., 2018). For a shallow and cool gas shale such as the Marcellus, gas adsorption in the matrix is important. Therefore, Eftekhari et al. (2020) modified the pressure diffusion equation to account for additional flow from adsorbed gas. Appendix Tables 10–12 show the key reservoir properties, gas composition, fluid properties, and adsorption parameters used by Eftekhari et al. (2020) to construct the master curve for the Marcellus Shale.where t˜=t/τ is the dimensionless time and c and a are the fitting parameters listed in the Appendix, Table 13. These parameters depend on gas and adsorption properties and drawdown pressures unique to each shale. Equation 4 was originally published as a three-parameter equation in equation B.1 in Saputra et al. (2021b). In the present paper, we successfully simplified it into a two-parameter equation by discovering that the product of the two fitting parameters a and b is exactly 1/2⁠. Consequently, if b=1⁠, t˜ is raised to the power of 1/2⁠, following the universal emergence of the square root of the time decline rate (Marder et al., 2018, 2021). Figure 9 plots each of the numerical master curves fitted with equation 4.where d is the half-distance between two hydraulic fractures, L is the half-length of hydraulic fractures (see Figure 8), and K is the initial slope of the master curve in Figure 9 or equation 4 with respect to t˜⁠, which is approximately 0.36 for the Marcellus. Because of the advancement of completion technologies in the Marcellus Shale, the values of KdL change over time as listed in the Appendix, Table 14.Figure 10 demonstrates the use of the physical scaling method to extend the historical well prototypes over up to two or three more decades. The thick black lines are the cumulative GEV means of annual gas rate (historical well prototypes) that were described in the previous section. After unit conversion from billion standard cubic feet to kilotons, we then match each of the historical well prototypes to the master curves in equation 4 or equation 5. The matching process is by scaling the x-axis, years on production, by the pressure interference time, τ⁠, and by scaling the y-axis, cumulative mass produced, by the mass of initial gas in place, M⁠. Table 3 lists the two matching parameters, τ and M⁠, for all 16 well cohorts.The physical scaling matches are shown in Figure 10 as the red and green lines. We determine EUR from each well cohort by looking at the endpoint of each physical scaling projection. We infer that the northeast core is the best region of the Marcellus, with EURs of up to 15 BSCF/well, followed by the southwest core with EURs of up to 11.5 BSCF/well. From a geological perspective (Dolson, 2016; Higley et al., 2019), the northeast core is more mature (dry gas) and has a thicker shale layer than the southwest core. Thus, its productivity is higher. The Marcellus noncore is less productive, with EUR lower by two- to threefold. The worst producing region is the Marcellus outer, which produces merely 1–2 BSCF/well, so that nobody has drilled in this area up until the present. From Figure 10, we conclude that the advancement of completion technologies helps to boost production almost twofold for both the northeast and southwest cores (comparing the oldest to newest cohorts). However, for the Marcellus noncore and Marcellus outer, where reservoir quality is poor, these technologies most likely do not help much.Well prototypes obtained from GEV statistics and physical scaling are highly idealized, with oil or gas production that lasts for decades. In reality, production from shale wells may not last that long in most cases. The reduction in reservoir pressure, increase in water cut, closure and deterioration of hydraulic fractures, and reduction of effective rock permeability decrease well production rates in shales.where Nactive,i and Ninactive,i are the numbers of active and inactive wells in year i. If we plot Psurvival for a well cohort in a region of a play versus year on production, and fit a parabolic curve, then the intercept of that curve at Psurvival = 0 defines the maximum time of well survival, tmax (months).Figure 11 displays the probabilities of well survival of the northeast core and southwest core wells (see SOM-2, supplementary material available as AAPG Datashare 173 at www.aapg.org/datashare, for all of the Marcellus regions). Each color represents a different year of completion. In both plots, the dotted line is the average of endpoint probabilities of survival of all of the well cohorts. The dashed line is a parabolic fit of the oldest well cohort that intercepts the maximum time of well survival, tmax, at which Psurvival = 0. We can see that wells in the northeast core survive longer than the wells in the southwest core. This may be due to different maturity levels of the reservoir rock. The southwest core is less mature, and we expect that there may be condensate banking near the wellbore that prohibits the wells from surviving for a longer time.Table 4 summarizes the average probability of well survival and maximum time of survival for each subregion of the Marcellus. We use these values to stop production of each well prototype shown in Figure 10.To confirm that we have obtained robust and stable well prototypes, we replace the production rate from each existing well with its corresponding well prototype, time shifted to the date of first production from the actual well, and sum up. The summation should match closely the historical gas rate. Figure 12 shows the result of this procedure for 12,406 existing wells in the Marcellus. The red and green curves are the physical scaling forecasts with and without exterior flow that match the black historical production curve. Both physical scaling projections forecast the play’s EUR of approximately 85 TSCF by 2030.We call this projection a base or do-nothing scenario of shale play development, in which there is no further drilling of new wells. This is the case for the already overdrilled, mature Barnett (Patzek et al., 2019) and Fayetteville plays, or if a global crisis hits the oil industry so hard that investment dries up and operators stop drilling new wells. In the Marcellus, if the operators stopped drilling today, then the total field gas rate would drop by half in just 3 yr. This is the true nature of shale play development, in which massive drilling must continue to maintain stable production. In the next three sections, we discuss in detail how we predict future infill drilling in the Marcellus.Figure 13 shows the distribution of tested cells in the northeast and southwest cores (see SOM-2, supplementary material available as AAPG Datashare 173 at www.aapg.org/datashare, for all of the Marcellus regions). The orange cells are locations where productive wells exist and the black cells are locations for dry hole wells.The number of productive and dry cells and the probability of success for each Marcellus subregion are summarized in Table 5. The highest probability of success exists in the northeast core of the Marcellus, whereas it is lowest in the Marcellus outer, as expected. However, these values are relatively low compared to a more mature shale play (e.g., Haynesville) (Saputra et al., 2021b), in which the probability of success is more than 90% (9 out of 10 wells drilled in Haynesville are expected to be productive). An immature shale tends to have a water saturation that is higher than those of oil or gas (Patzek et al., 2019; Saputra et al., 2019). If the effective permeability to water is higher than those to hydrocarbons, then the wells are expected to flow little or no oil or gas, leading to dry wells.In this section, we calculate the number of potential infill wells to be drilled. This calculation is similar to the previous one. First, we cover each subregion of the Marcellus with the 1-mi2 grid cells. Due to regulatory restrictions, it is impractical to drill and produce inside-city boundaries. Therefore, we exclude every grid cell that intersects urban areas, shown as the dark gray areas in Figure 3.where Ngrids is the total number of 1-mi2 grids in each subarea, ψ is the difference between the maximum allowable number of wells per square mile to avoid fracturing hits and the number of existing wells per grid cell, and Psuccess is the probability of success that we calculated in the previous section.Figures 14 and 15 help demonstrate the calculation of potential infill wells in the Marcellus by showing the lateral directions, well density, and infill potential maps of the northeast and southwest cores of the Marcellus (see SOM-2, supplementary material available as AAPG Datashare 173 at www.aapg.org/datashare, for all subregions of the Marcellus). From Figures 14A and 15A, we observe how operators adapt to the direction of minimum horizontal stress in the Appalachian Basin. To generate stable transverse vertical hydraulic fractures, the operators should drill wells at ∼S20°E in the northeast core and ∼S45°E in the southwest core.For each grid cell, we count how many wells intersect this cell in the subsurface. The outcome is the subsurface well density maps shown in Figures 14B and 15B. The different colors show how many wells already exist per square mile. Following Saputra et al. (2019, 2021b), we assume that the maximum allowable number of wells to be drilled is 4 wells/mi2 for a typical fracture length of 1200 ft. Based on equation 8, we calculate ψ by subtracting each well density map from the assumed maximum allowable number of wells. Figures 14C and 15C show the results as (|ψ|+ψ)/2⁠. Notice that a negative ψ converts to zero in both maps, suggesting overdrilling in each black grid cell.Finally, we sum the potential infill wells for all grid cells and multiply the result by the probability of success. The results are tabulated in Table 6.Based on the infill potentials in Table 6 and historical drilling rates in Figure 16A, we schedule the future drilling programs shown in Figure 16B. We have obtained the number of active rigs (rig count) in the Appalachian Basin from EIA, and the number of horizontal wells completed per month (drilling rates) from Enverus (2021). It is interesting to see that after 2020, the monthly rig counts match perfectly the drilling rates. We also highlight three levels of drilling rates: 130 wells/month in 2011 to 2015, with gas prices above $10/MSCF; 80 wells/month in 2017 to 2019, with gas prices recovering after the 2008–2009 financial crisis to >$3/MSCF; and 50 wells/month during the oil crisis of 2016, followed by the global coronavirus disease 2019 pandemic of 2020 to 2021.Based on Table 2, we observed that operators mostly drill in the core areas. Therefore, we continue this trend and schedule to drill first the 2406 and 1278 potential wells in the northeast and southwest cores, respectively, using the historically low drilling rate of 50 wells/month. With this infill drilling program in place, in 2027, it is expected that no well locations will be left in the core areas. Thus, operators may be forced to move to the less productive Marcellus noncore and drill all of the potential 7896 wells at a higher drilling rate of 80 wells/month. Here, we assume that the industry will recover by 2027. By 2035, there should be no new locations left for infill drilling in the core and noncore areas of the Marcellus. Thus, operators may be able to drill only the least productive Marcellus outer area at the most optimistic drilling rate of 130 wells/month. In the next section, we demonstrate that drilling ∼15,000 wells in the unproductive outer area is economically unfeasible, except at high gas prices.Using an approach similar to that in the previous section, we replace the future wells scheduled in the drilling programs in Figure 16B with their corresponding time shifted well prototypes and sum up the resulting rates to obtain the total field rate. In Figure 17A, we show the total field rate as the increments from the base or do-nothing case in Figure 12. We gradually infill both core areas and then the noncore area; finally, we infill the Marcellus outer area. By drilling in the core areas, we can maintain a production plateau of approximately 25 BSCF/day for the next 5 yr. Even if we increase the drilling rate from 50 to 80 wells/month, the infill in the noncore area does not prevent production from declining. The most optimistic and highly unlikely infill scenario of 130 wells/month in the outer area still cannot save the Marcellus from its terminal decline.Finally, we integrate the total field rates in Figure 17A to predict the ultimate recoveries from the Marcellus (see Figure 17B). The 12,406 existing wells in the Marcellus are predicted to ultimately produce 85 TSCF of gas by 2040. By adding the 3684 potential wells in the core areas, the cumulative production is poised to increase to ∼150 TSCF by 2040. The additional 7896 wells in the noncore area may increase the ultimate production to only ∼180 TSCF by 2050. Lastly, the staggering number of 15,000 wells in the outer area may yield only an additional ∼20 TSCF, an endeavor that is economically unfeasible, as we discuss in the next section.In Table 7, we summarize the existing production, proven reserves, remaining resources, and EUR from four different reservoir qualities of the Marcellus. Besides the expected values (mean), we also present the upper bounds (⁠P10⁠) and the lower bounds (⁠P90⁠) of each estimated reserve–resource based on GEV statistics. We also compare the results of this study with other published results (Ikonnikova et al., 2018; Higley et al., 2019; Energy Information Administration, 2020a) (see Table 8). Our estimate of the remaining resources at 88 to 109 TSCF is highly similar to what USGS predicted in 2019, at 96.5 TSCF. Next, we translate the value of the remaining resources into EUR or TRR by adding existing production and proven reserves. Our predicted TRR value at 178 to 200 TSCF is much lower than what EIA and BEG predicted at 310.6 and 560 TSCF, respectively.To clarify, our analysis was completed solely for the Marcellus. The USGS prediction is also based on the Marcellus layer at a similar target depth, and this is probably why our estimates are so similar. There may be additional potential production coming from the adjacent formations, such as the Upper Devonian Marcellus. (However, as of January 2021, there were only 399 wells producing from the Upper Devonian interval, and their production was negligible compared with the ≥15,000 wells drilled in the Marcellus.) Perhaps EIA and BEG also considered these additional layers, and therefore their estimates were two to three times higher than ours. The lower ultimate resource estimates proposed in this paper are supportable, we think, from the high statistical confidence we have in well productivity obtained with the approach taken.where tmax is the lifetime of the project. In this case, it is the maximum time on production that we have calculated from the probability of well survival in a given area of the Marcellus. A profitable project should yield a positive NPV, whereas a negative NPV means a money loss. When NPV=0⁠, we call it breakeven, which occurs at a certain gas price. Our goal is to calculate the NPV of each drilling project and determine the gas price at which the breakeven of that project happens. The lower the breakeven price, the better the project.Table 15 in the Appendix lists the economic parameters used to calculate the NPVs of different development scenarios in the Marcellus Shale. We have obtained these parameters from various sources (Hefley and Seydor, 2011, 2015; Duman, 2012; Khodabakhshnejad et al., 2019; Range Resources, 2019; Tabuchi, 2020). The gross revenue is the summation of annual gas and natural gas liquid (NGL) production multiplied by the gas and NGL prices, respectively. Using GEV statistics, we calculate the expected value of condensate to gas ratio that can be used to estimate the amount of NGL produced (see the Appendix, Table 16). The capital expenditure is the summation of drilling and completion cost (∼$5 million), land acquisition cost (∼$0.5 million), and plug and abandonment cost (∼$0.3 million). The operating expenditure is the total annual gas and NGL production (in barrels of oil equivalent [BOE]) multiplied by the operating cost (∼$4.2/BOE). The royalty is the gross revenue multiplied by the royalty rate (∼15%/yr). The overall tax rate is the most difficult parameter to calculate (see SOM-3, supplementary material available as AAPG Datashare 173 at www.aapg.org/datashare for details). To calculate taxes, we need to assume the severance tax rate (∼$0.3/BOE), corporate tax rate (∼25%/yr), and intangible expenditure (∼50%/yr). Finally, the 10% DIS is chosen because it is widely used in the annual reports of United States oil and gas companies.In the present study, we use two scenarios of NPV calculations for four reservoir qualities in the Marcellus. We show the results in Figure 18. In the first scenario, we fix the NGL price at $31.5/bbl as of August 2021 (Cocklin, 2021b) while varying the gas price to determine at what gas prices the infill project breaks even. In the second scenario, we fix the gas price at $2.81/thousand SCF (KSCF) as of August 2021 (Cocklin, 2021a) and vary the NGL prices. The detailed NPV calculations at different gas and NGL prices are tabulated in SOM-3 (supplementary material available as AAPG Datashare 173 at www.aapg.org/datashare).In Figure 18A, we see that both core areas (northeast and southwest) have similar breakeven points at a gas price less than $1.8/KSCF. The Marcellus noncore is less productive than the core areas, and thus the breakeven point shifts to $3.1/KSCF. As expected, the least productive Marcellus outer has the highest breakeven point at $7/KSCF. For comparison, a comprehensive study by the Journal of Petroleum Technology (2019) shows similar breakeven prices for the Marcellus, ranging from $1 to $8/KSCF.In Figure 18B, we see the impact of NGL prices on the profitability of each infill project. The northeast core contains only dry gas with no NGL production, so that the NPV remains constant whatever the NGL price. However, the southwest core yields a significant amount of NGL, so that the higher the NGL price, the more profitable the infill project. Lastly, both Marcellus outer noncore and Marcellus outer give negative NPV for all NGL prices shown in the graph, and they are not profitable at current prices.From this exercise, one can infer that drilling in both the Marcellus northeast and southwest core areas may always yield a positive cash flow, even if gas prices plunge below $2/KSCF. The Marcellus noncore area may be an alternative for infill drilling if gas prices exceed $3/KSCF. Finally, we do not recommend drilling in the Marcellus outer area because it is expected to generate negative cash flows, unless gas prices increase to more than $7/KSCF, which is unlikely in this era of low oil and gas prices.We have provided an optimal play-wide assessment of the Marcellus Shale by considering the shale play geology, advancement of well completion technologies, physics of natural gas production from the horizontal hydraulically fractured wells, and economics of drilling projects.Our GEV statistics approach makes all of the well prototypes robust and in excellent agreement with the physics-based scaling curves. Using these prototypes, we were able to match rather well the historical gas production from the entire Marcellus Shale.Both core areas in the northeast and the southwest of the Marcellus give the highest EUR and the most profitable wells due to the more mature and thicker Marcellus Formation.The newer wells yield higher EURs due to (1) longer lateral lengths, (2) larger hydraulic fractures, and (3) more hydraulic fracture stages. However, in areas of poor reservoir quality, these advancements in completion technologies cannot help much.Ultimately, the 12,406 existing wells in the Marcellus are expected to produce 85 TSCF of natural gas by 2040. By adding another 3864 and 7896 potential wells in the core and noncore areas, the ultimate recovery is poised to increase to 150 and 180 TSCF, respectively.Operators should not bother with drilling in the outer area of the Marcellus because it may be unprofitable for all of the scenarios we considered. Other attractive shale plays should be drilled that may yield higher cash flows (e.g., Haynesville Formation and Utica Shale).In our opinion, the hybrid data-driven and physics-based approaches are the future of production forecasting and reserve estimation in all shale plays. These methods deliver the best—in the least squares sense—predictions of possible futures, are free from bias, and avert significant over- or underpredictions of reserves.The authors thank KAUST for supporting this work through baseline research funding to Tadeusz Patzek. The Computer, Electrical, and Mathematical Sciences & Engineering Division at KAUST supported Wissem Kirati. We also thank the AAPG Bulletin technical editor, Andrea Sharrer, and all four reviewers for their thorough, informative, and constructive reviews that have helped us to improve substantially the manuscript. In particular, we thank the editor and reviewer 4, who excelled in their detailed, insightful remarks that were a pleasure to implement.Supporting online materials 1–3 are available in an electronic version on the AAPG website (www.aapg.org/datashare) as Datashare 173.Wardana Saputra is a research associate at the Hildebrand Department of Petroleum and Geosystems Engineering, The University of Texas at Austin. He received a B.Sc. degree in petroleum engineering from Institut Teknologi Bandung, Bandung, Indonesia, in 2015. From the King Abdullah University of Science and Technology in Thuwal, Saudi Arabia, he earned an M.Sc. degree in earth science and engineering and a Ph.D. in energy resources and petroleum engineering in 2017 and 2021, respectively. During his Ph.D., he analyzed a huge data set of ≥500,000 shale–tight oil and gas wells in all of the major petroleum basins in the United States to create a new physics-based data-driven production forecasting method in shales.Wissem Kirati was a research engineer at the Ali I. Al-Naimi Petroleum Engineering Research Center at KAUST, Thuwal, Saudi Arabia. He has several years of experience working in the oil and gas field worldwide. He was also involved in consultancy-based studies for major operating oil companies in drilling, reservoir, and research and development areas. Before joining KAUST, he was a drilling engineer at Schlumberger and a reservoir engineer at Euro Engineering AG.David Hughes is an earth scientist who has studied the energy resources of Canada and the United States for more than four decades, including 32 years with the Geological Survey of Canada as a scientist and research manager. He is the president of Global Sustainability Research, a consultancy that has analyzed the geological fundamentals and production potential of unconventional oil and gas plays across Canada and the United States. He has published and lectured widely on energy and sustainability issues in North America and internationally. He is also a fellow of the Post Carbon Institute, a board member of the Physicians, Scientists, and Engineers for Healthy Energy, and a research associate with the Canadian Centre for Policy Alternatives.Tadeusz (Tad) W. Patzek is the energy resources and petroleum engineering professor and the center director of the Ali I. Al-Naimi Petroleum Engineering Research Center at KAUST. Before coming to KAUST, he was the Lois K. and Richard D. Folger Leadership Professor and chairman of the Petroleum and Geosystems Engineering Department at The University of Texas at Austin. He also held the Cockrell Regents Chair #11. Between 1990 and 2008, he was a professor of geoengineering at the University of California, Berkeley. Before joining Berkeley, he was a researcher at Shell Development. Patzek’s research involves mathematical and numerical modeling of earth systems with an emphasis on subsurface fluid flow. He is working on the thermodynamics and ecology of human survival and energy supply for humanity. He is also working on unconventional natural gas and oil resources. Patzek is a coauthor of more than 400 papers and reports and 1 book, and is currently writing 2 books. In 2020, he received the highest technical honor bestowed by the European Association of Geoscientists and Engineers, the Desiderius Erasmus Award for lifetime contributions. In 2022, Patzek was awarded the Society of Petroleum Engineers IOR Pioneer Award.