Generalized Extreme Value Statistics, Physical Scaling and Forecasts of Oil Production in the Bakken Shale

Abstract

We aim to replace the current industry-standard empirical forecasts of oil production from hydrofractured horizontal wells in shales with a statistically and physically robust, accurate and precise method of matching historic well performance and predicting well production for up to two more decades. Our Bakken oil forecasting method extends the previous work on predicting fieldwide gas production in the Barnett shale and merges it with our new scaling of oil production in the Bakken. We first divide the existing 14,678 horizontal oil wells in the Bakken into 12 static samples in which reservoir quality and completion technologies are similar. For each sample, we use a purely data-driven non-parametric approach to arrive at an appropriate generalized extreme value (GEV) distribution of oil production from that sample’s dynamic well cohorts with at least $1,2,3,\cdots$ years on production. From these well cohorts, we stitch together the $P_{50}$, $P_{10}$, and $P_{90}$ statistical well prototypes for each sample. These statistical well prototypes are conditioned by well attrition, hydrofracture deterioration, pressure interference, well interference, progress in technology, and so forth. So far, there has been no physical scaling. Now we fit the parameters of our physical scaling model to the statistical well prototypes, and obtain a smooth extrapolation of oil production that is mechanistic, and not just a decline curve. At late times, we add radial inflow from the outside. By calculating the number of potential wells per square mile of each Bakken region (core and noncore), and scheduling future drilling programs, we stack up the extended well prototypes to obtain the plausible forecasts of oil production in the Bakken. We predict that Bakken will ultimately produce 5 billion barrels of oil from the existing wells, with the possible addition of 2 and 6 billion barrels from core and noncore areas, respectively.

1. Introduction

Over the last decade, crude oil and gas production from hydrofractured shales in the United States has accounted for most of the net increase of global crude oil production. Therefore, it is important to have a reliable, quantitative method for delineating the possible futures oil and gas production in the data-rich US shales. The current industry-standard methods of forecasting production from shales are variants of the empirical decline curve analysis (DCA) [1], developed 75 years ago. Lately, some of the more sophisticated methods, for example, Fractional Decline Curve [2] have become popular.

Unlike other analytical and numerical methods that require numerous reservoir parameters and a lengthy calculation or simulation time, DCA only requires production data to predict future production by extrapolating oil or gas rate observed in a boundary-dominated flow regime. Because of DCA’s simplicity, most petroleum engineers adopt it for reserve assessment in shales, including Estimated Ultimate Recovery (EUR) predictions from USGS [3,4] and EIA [5].

USGS first split the Bakken region into 6 assessment units (AUs) and defined sweet spots. Then, they calculated the number of wells that could be drilled in each AU by dividing its total area with the average drainage area per well. In parallel, they used DCA to calculate the average EURs for sweet spots and other areas. The 7.4 billion barrels of undiscovered technically recoverable oil was obtained by multiplying the total number of wells that could be drilled in each AU with the corresponding average EUR multiplied by the untested fraction of that AU. EIA used a similar approach by dividing the Bakken region into 5 AUs and refining them by counties to determine the infill well potential. Both USGS and EIA predictions were assessed by Hughes [6,7,8,9]. We predict that the undiscovered, technically recoverable oil in the Bakken is 2 billion bbls in the core area and it might be 6 billion bbls in the noncore area (Figure 1) or 8 billion barrels in total.

Most shale wells do not reach the boundary-dominated flow regime for their entire production lives because of the vanishing matrix permeability. Thus, the traditional DCA frequently overestimates EURs of shale wells. To address this issue, many authors have suggested improved DCA methods, specific to shale wells: the Power Law Exponential Decline [10], the Stretched Exponential Decline [11], the Logistic Growth Model [12], the Extended Exponential DCA [13], and the Extended Hyperbolic DCA [14]. To make things worse, the empirical DCA fits of particular wells are ill-suited to forecasting production from a wide area of a given shale play in which reservoir properties vary and uncertainties abound. Therefore, some authors have developed probabilistic models to introduce a range of possible outcomes into their production forecasts [15,16,17,18,19,20]. The most common assumption is that well productivities in shales are log-normally distributed.

In this paper, we adopt a hybrid data-driven and physics-based method of predicting oil or gas production in shales that has been introduced in our previous work [21,22]. Here, we consider only black oil production. First, we identify play regions in which reservoir quality is similar, see Figure 1. In each region, we identify well classes by different completion technologies. Finally, a well class in a region constitutes a well sample. We ensure that oil production from all wells in each sample is statistically uniform, that is, has a unimodal distribution. For each well sample, we then identify well cohorts with at least $1,2,\cdots$ years on production. In general, well cohorts contain different sets of wells that satisfy the minimum time on production required for each cohort. It turns out that each cohort of wells is superbly characterized by its unique Generalized Extreme Value (GEV) distribution (see Appendix A) of annualized well rates or cumulative well production. Different cohorts in the same sample have different GEV distributions, each with its unique expected value, median and mode. Here we choose the somewhat better GEV fits of the production rate distributions. Each GEV distribution is statistically superior to the corresponding log-normal distribution at the 95% confidence level. When we plot the expected values of the GEV distributions of all wells cohorts in a sample versus elapsed time of production, we obtain this sample’s average $P_{50}$ statistical well prototype that is purely field data-driven.

Now we fit each statistical well prototype with a physical scaling curve that extends this prototype to 30 years on production. The physical scaling curves are based on an analytical solution of the pressure diffusion equation in the hydrofractured horizontal shale well geometry. In previous papers, we comprehensively detailed the physical scaling solutions for shale gas wells [23,24] and shale oil wells [22,25]. Late-time flow from outer reservoir encompassing the stimulated reservoir volume (SRV) [26] was also quantified. We have verified that our physical scaling is equivalent to detailed numerical reservoir simulations [27]. This scaling is much simpler to set up and runs almost infinitely faster than the corresponding reservoir simulations.

At this point, for each well sample in every play region, we have obtained a unique hybrid mean ($P_{50}$) well prototype with 30 years on production. In each play region, we know how many wells were drilled and completed each month up until current time. We then multiply each well prototype by the number of wells completed per month and stack them up to represent the total historical field rate and future production decline. In this manner, we obtain a ‘base’ case forecast for all existing wells. This base case forecast is a ‘do nothing’ scenario with no new wells drilled in the future. For all other forecasts of future field production rate, we first determine the infill potential or the number of wells that can be drilled in the future without causing significant interwell pressure interference (fracture hits). We cover each region with a fish-net grid that consists of one-square-mile pixels. We then calculate the infill potential as the number of wells that can be drilled so that the total number of wells in each pixel is less than the maximum allowable number of wells without fracture hits. Next, based on the infill potentials for all regions, we create future drilling programs to obtain plausible forecasts of oil or gas production. Based on current rig count in the Bakken, we assume a constant overall drilling rate. Finally, we assign the correct well prototype to every future well that will be drilled during each month of a postulated drilling schedule, and sum them up to obtain a forecast scenario.

In this paper, we select the Bakken shale, the current second-largest oil producer in the U.S. with 1.5 million bbl of oil per day, as an illustration. Being one of the oldest shale oil plays, Bakken has been a field laboratory to test drilling and completion technologies and increase well productivity. Currently, Bakken has ∼15,000 active hydrofractured horizontal wells with a few wells that have 18 years of production data. In a previous paper [22], we scaled well-by-well all ∼15,000 wells in the Bakken. We accounted for well refracturing and/or changes in downhole pressure. It turns out that the 12 well prototypes obtained with our hybrid GEV–physical scaling method are as good in duplicating the total field rate as the super-precise scaling of each individual well in our previous work [22]. Given the results of our analyses that are free of bias, policy-makers should not assume that the production boom in the Bakken shale will last decades longer.

2. Results

Our approach is as follows:

• Divide all 14,678 horizontal oil wells in the Bakken shale into 12 samples in which oil production is statistically uniform;
• Fit a generalized extreme value distribution to all wells in every sample and obtain 12 stable mean $P_{50}$ well prototypes;
• Fit the physics-based scaling curves to every statistical well prototype and extend these prototypes smoothly to 30 years on production;
• Replace oil production from all existing wells with the 12 extended well prototypes and obtain a ‘base’ forecast;
• Calculate the infill potential for each of the 12 regions in the Bakken; and
• Create the plausible infill drilling schedules and forecast total field oil production rate up to the year 2050.

2.1. Design of Well Regions

In the previous work [21], we divided all 13,141 horizontal wells in the Barnett shale into 6 static regions that corresponded to the geographical borders of 6 major gas-producing counties. We justified this approach as follows: (i) Productive Barnett shale is a single gas reservoir; (ii) the Barnett play area is relatively small, and the division into counties captures different reservoir qualities; and (iii) completion technologies have not been revolutionized over time.

However, in the Bakken shale, we introduce a more complicated set of 12 static well samples classified by reservoir, geography, and completion dates. In Bakken, there are two different producing reservoirs, the Middle Bakken and Three Forks. Accordingly, we first split all Bakken wells into the Middle Bakken and Three Forks groups. To account for different reservoir quality, we then divide each of these two well groups into two macro regions: the core and noncore areas, see Figure 2. Finally, we further refine the four Bakken macro regions by splitting them into three classes by completion date intervals, 2000–2012, 2013–2016 and 2017–2019. These classes reflect advancements in well completion technology, such as longer wells, more fracture stages, fewer clusters and perforation shots, bigger hydrofractures, more proppant, and so forth. In the end, we divide the existing 14,678 horizontal oil wells in the Bakken into 12 unique samples listed in

Table 1. Static regions for the Bakken shale play.

Well	Reservoir	Area	Completion	Sample
Sample			Date	Size
1	Middle Bakken	Noncore	2000–2012	2550
2	Middle Bakken	Noncore	2013–2016	1355
3	Middle Bakken	Noncore	2017–2019	223
4	Three Forks	Noncore	2000–2012	735
5	Three Forks	Noncore	2013–2016	1204
6	Three Forks	Noncore	2017–2019	207
7	Middle Bakken	Core	2000–2012	2086
8	Middle Bakken	Core	2013–2016	2534
9	Middle Bakken	Core	2017–2019	1112
10	Three Forks	Core	2000–2012	428
11	Three Forks	Core	2013–2016	1502
12	Three Forks	Core	2017–2019	742
TOTAL				14,678

.

2.2. Statistical Well Prototypes

For each of the 12 Bakken well samples defined in Section 2.1, we construct a statistical well prototype based on the Generalized Extreme Value (GEV, see Appendix A) mean or the $P_{50}$ values of annual rates from all well cohorts that exist in that particular region. Figure 3 illustrates the procedure of arriving at the GEV well prototypes. For every GEV distribution fit we acquire values of the location parameter ($\mu$), scale parameter ($\sigma$), shape parameter ($\xi$), confidence interval (CI), mean ($P_{50}$), median, mode, $P_{10}$ and $P_{90}$. In this paper, we show only two GEV fit examples (Figure 4 and Figure 5). Please see Supporting Online Materials-1 to find all GEV fits for all 12 regions of the Bakken shale. The resulting statistical well prototypes are plotted versus the square root of time on production in Figure 6. In the initial year on production, for some well cohorts, time was shifted by subtracting the first 1–3 months on production (see the second column of

Table 2. Scaling parameters for an average well in the Bakken regions.

Well Sample	Physical Scaling				Physical Scaling with Exterior Flow
	Shift Years	$\mathit{\tau }$ Years	$\mathcal{M}$ kbbl/Well	EUR kbbl/Well	$\mathit{\tau }$ Years	$\mathcal{M}$	a	EUR kbbl/Well
Middle Bakken noncore [2000–2012]	0.25	24.5	454.3	271.1	19.9	286.7	0.050	282.2
Middle Bakken noncore [2013–2016]	0.00	11.8	341.3	211.7	10.9	267.2	0.026	236.4
Middle Bakken noncore [2017–2019]	0.00	5.0	346.9	215.5	5.0	306.6	0.015	259.5
Three Forks noncore [2000–2012]	0.15	13.3	288.4	207.5	10.7	181.2	0.058	238.3
Three Forks noncore [2013–2016]	0.00	10.8	271.9	196.1	9.9	212.4	0.030	221.1
Three Forks noncore [2017–2019]	0.00	5.0	270.8	195.4	5.0	238.4	0.018	236.5
Middle Bakken core [2000–2012]	0.20	25.0	810.8	482.5	22.6	536.3	0.050	512.2
Middle Bakken core [2013–2016]	0.00	9.1	543.8	337.7	8.2	420.7	0.026	389.7
Middle Bakken core [2017–2019]	0.00	5.0	581.9	361.5	5.0	513.6	0.015	434.6
Three Forks core [2000–2012]	0.27	25.0	681.1	470.9	25.0	473.4	0.058	511.6
Three Forks core [2013–2016]	0.00	10.3	450.1	324.6	9.6	353.1	0.030	369.7
Three Forks core [2017–2019]	0.00	5.0	462.8	334.0	5.0	406.9	0.018	403.5

) to discount all possible initial disturbances in well performance. After these time shifts, each mean or $P_{50}$ well prototype starts from a straight line versus the square root of time on production, as expected in linear flow regime. Later, cumulative oil production bends down due to the inter-hydrofracture pressure interference and exponential decline of production rate [23]. For our GEV distribution fits, the mean or $P_{50}$ prototype is always higher than the median, and the median is higher than the mode. The $P_{10}$ and $P_{90}$ prototypes diverge from the mean with time, indicating that with time the best wells grow better and worst wells worse. A similar trend was also observed in the Barnett shale [21].

2.3. Physical Scaling Fits

For each region, we fit the statistical $P_{50}$ well prototype with a physical scaling curve, see Appendix B, to extend this prototype to 30 years on production. The fitting procedure is detailed in the Materials and Methods section. Briefly, by least squares fit, we find optimum values of the scaling parameters $\tau$ and $\mathcal{M}$ so that the physics-based well prototypes match the statistical well prototypes. The results are shown in Figure 6. The red solid lines show the most pessimistic versions of each physics-based prototype, which assume that oil is produced only from the interior of the stimulated reservoir volume (SRV). The green solid lines show a more realistic forecast with an assumption that at late times there will be additional exterior flow of $a\sqrt{t}$ towards the SRV. Finally, one obtains the Estimated Ultimate Recovery (EUR) from each region by identifying the endpoint of cumulative production for the physics-scaling prototype. In Table 2, the scaling parameters, $\tau$, $\mathcal{M}$ and a, and EUR values are listed for each region with and without exterior flow. Notice that in general the newest wells have the shortest pressure interference times, $\tau$, and fastest decline rates.

2.4. Base or ‘Do Nothing’ Forecasts

We have replaced actual field production rate from all existing wells in each region with this region’s well prototype multiplied by the historic numbers of completed wells, and we time-shifted the products. The superposed prototypes match the historical field rate. Figure 7 shows the sum of the well prototypes with and without exterior flows (green and red lines) versus the total field historical oil rate and cumulative oil (black lines). The results are rather satisfactory. Our GEV prototypes are robust. They capture all peaks of total historical oil rate and flawlessly match the cumulative production. As stated before, the red line forecast is pessimistic. It assumes that there will be no production contribution from the reservoir exterior to the SRVs of individual wells. This pessimistic case may hold if a shale play, here Bakken, is already overdrilled. For the Bakken shale, we can infer that the 14,678 existing wells will ultimately produce 4.5–5.3 billion of oil by 2050. This is a ‘base’ or ‘do nothing’ case with no future drilling in the Bakken shale.

2.5. Infill Potentials

We have calculated the number of wells that can be drilled in the future in every one-square-mile pixel of the grid that covers the entire Bakken play. In order to calculate infill potentials, one should first determine well density. However, the publicly available data rarely provide information about the bottomhole locations of the wells. Instead, only surface locations are reported as latitude-longitude coordinates. Therefore, we have developed an algorithm that allows us to predict the bottomhole well density from surface well locations. See Figure 8 and the Materials and Methods section for more details. As mentioned before in Section 2.1, the Bakken shale has two reservoirs, the Middle Bakken and Three Forks and two significantly different reservoir qualities, the core and noncore areas. Thus the infill potentials should be calculated for these four parts of the Bakken. Supporting Online Materials-2 shows the gridding, well locations, wellhead densities, well densities and infill potentials for the four Bakken parts.

Table 3. Infill potential for Bakken region.

Reservoir	Region	Total Area	Existing	Infill
	Type	sq Miles	Wells	Potential
Middle Bakken	Noncore	8534.5	4128	14,171 *
Three Forks	Noncore	6905.9	2146	12,341 *
Middle Bakken	Core	3382.5	5732	2577
Three Forks	Core	1914.1	2672	1825
TOTAL		20,737	14678	30,914

* Given the low reservoir quality in noncore areas, our infill potentials there are probably too high.

lists the numbers of wells that can potentially be drilled in the future. We note in passing that infilling the less productive noncore areas to the same well density as the core areas is likely too optimistic.

2.6. Future Drilling Scenarios and Infill Forecasts

We have created several future drilling scenarios based on infill potentials listed in Table 3. We assume that the rig count in the Bakken will not significantly change from the current value and 120 wells will be drilled each month between now and 2041 (Figure 9a). By looking at the current trend, many operators have narrowed their drilling choices to only the core areas to avoid spending money on the less productive wells with high watercuts typical of the noncore areas [22]. Thus, it is reasonable to create a drilling schedule that exhausts all potential wells in the core areas first, before moving out to the less productive noncore areas. Figure 9b shows this scenario predicting that the Bakken core area will be completely drilled out by 2022. Later, the drilling in noncore areas can last up to 2041, according to our calculation. See Supporting Online Materials-2 for the more detailed schedules.

Next, we assign each new well scheduled to be drilled in the future to its corresponding region and class prototype. The results are shown in Figure 10. The black solid lines are the historical total field oil rate and cumulative oil. The purple lines denote the ‘base’ forecast from all 14,678 existing wells displayed in Section 2.4. The red lines show the result of adding 4402 new wells in the core area. Because all wells in the core area are high-grade wells, the total oil production rate climbs up its highest production peak at 1.6 million bbl/d in the year 2022. After 2022, there is no further drilling in the core area and the production rate declines by a factor of two within one year. Ultimately, the total of 19,080 wells (existing + new core) will produce 7 billion bbl of oil by the year 2050. The blue lines show the most optimistic case with not only the core areas but also the noncore areas will be drilled out by the operators. The drilling of less productive wells in the noncore areas will maintain the plateau oil rate of 1 million bbl/d up to the year 2041. However, this plateau will require the drilling of additional 26,500 wells (almost twice the current number of wells) which will have high water cuts. Ultimately, the total of 45,580 wells (existing + core + noncore) will produce 13 billion bbl of oil by 2050.

3. Discussion

We have presented an alternative to the current industry-standard empirical forecasts of oil production from hydrofractured horizontal wells in shales. With our hybrid modeling approach, we have matched current oil production in the Bakken rather accurately. We have also delivered an optimal prediction of possible futures of the Bakken shale play for up to three decades.

Our Bakken oil forecasting method extends the previous work on predicting fieldwide gas production in the Barnett shale [21] and merges it with our new scaling of oil production in the Bakken [22]. Our field data-driven statistical well prototypes are conditioned by well attrition, hydrofracture deterioration, pressure interference, well interference, progress in technology, and so forth. With no physical scaling, these prototypes follow the exact physics of linear transient oil flow with pressure interference. Therefore these statistical well prototypes serve as templates to calibrate the parameters of our physical scaling model ($\tau$ and $\mathcal{M})$ [23,24] and obtain a smooth time-extrapolation of oil production that is mechanistic and not merely an empirical decline curve. At late times, we add to our extended prototypes some radial inflow from the outside of well SRVs [26].

The extended $P_{50}$ well prototypes in Figure 6 can be used to compare ultimate recovery in each of the static regions we have identified in the Bakken. The lower bounds are the extended $P_{50}$ prototypes without exterior flow (red lines). In most cases, wells completed in the upper Three Forks reservoir are somewhat less productive than those in the Middle Bakken reservoir. The reasons for this difference are: (1) higher water saturation and water cut in Three Forks, (2) faster decline rate (lower $\tau$) and (3) lower initial oil in place (lower $\mathcal{M}$) [22]. This difference is consistent with the stratigraphic column of the Bakken total petroleum system, where Middle Bakken is sandwiched between two world-class source rocks with 10% TOC, the Upper Bakken Shale and the Lower Bakken Shale. On the other hand, the Three Forks formation is below the Lower Bakken Shale member and is exposed to water-bearing formations beneath [3,28,29,30]. For the same reasons, the core and effective drilling areas in the Three Forks are smaller than those in the Middle Bakken.

In both reservoirs, production from the core areas is superior to that from the noncore areas. The core area located in the center of Williston Basin has been known as the most oil-prolific location in the Bakken region in North Dakota [3,28,29,30]. Since the 1950s, oil has been produced there in the thickest, naturally fractured Middle Bakken formation in the Nesson anticline. One inexpensive vertical well drilled in the core area in the 1950s has had ultimate recovery of 200 kbbl, the same as a $10 million hydrofractured horizontal well drilled in the Middle Bakken noncore area nowadays. The noncore areas are less productive because the three Bakken formations: Upper, Middle and Lower Bakken are pinching out upward (are thinner and less mature) near the edges of the Williston basin. Consequently, the noncore areas are producing more water than oil, with watercut exceeding 50% on average.

The newly completed wells have much higher initial oil rates than the older ones, because they have: (1) longer lateral lengths, (2) bigger hydraulic fractures and (3) more fracture stages [22,31,32,33,34]. However, the newly completed wells decline faster and have essentially the same ultimate recovery as the older wells. The reasons for this behavior have been described elsewhere [35]. Interestingly, in most cases, older wells completed in 2000–2012 have higher ultimate recovery than newer ones, even though their initial production rates are lower. These older wells might have been drilled in the best locations ever in the Bakken region. In addition, shorter lateral lengths and fewer fracture stages may help in maintaining a stable pressure drawdown and prevent reservoir degassing that is unfavorable for future production. For comparison, the average lateral length in 2005 was 5000 ft while the average lateral length in 2019 doubled to 10,000 ft. Historically, the number of hydrofracture stages in the Bakken has increased over time from 8 stages in 2007 [36] to 18 stages in 2009 [31], 35 stages in 2016 [32] and to as many as 60 stages in 2019 [33].

According to our records, more than 90% of the wells completed after 2017 are located in the core areas only. Operators have learned to drill only the best parts of the Williston Basin and avoid the less mature noncore areas. However, after calculating the infill potentials of all areas, we predict that by 2021 there will be no well locations left for future drilling in the core areas. Assuming a constant current drilling rate of 120 wells per month, the total field oil rate in the Bakken will reach record level of about 1.6 million bbl/d in 2021. Without further drilling, production will decline by one-half within a year. Later, operators will be forced to drill in the less productive, high watercut noncore areas along the edges of the Williston Basin. Our findings suggest that policy-makers should not assume that the shale oil boom in the Bakken will last for several decades longer. We recommend that operators not focus only on increasing the initial oil rate. Maintenance of reservoir pressure above the bubble point by preventing over-drilling is key to increasing ultimate oil recovery.

4. Materials and Methods

We mined public well data in the Bakken region from DrillingInfo database, the Montana Department of Natural Resources and Conservation–Board of Oil and Gas Conservation website and the North Dakota Department of Mineral Resources–Oil and Gas Division website. From DrillingInfo, we selected only 14,678 horizontal oil wells that were completed in the Middle Bakken and Three Forks formations after January 2000. We designed an integrated MATLAB${}^{\circledR }$ software package to perform data cleanup, consolidation and unit conversions. The clean production data for each well consist of a vector of elapsed times on production, vector of oil and water production rates and records of wellhead location (latitude and longitude), lateral length, completion date, reservoir description and maximum oil rate. The average reservoir and well lateral properties are listed in Appendix C. Our approach can be divided into six main steps:

• Define play regions in which oil production is statistically uniform. In each region, follow the steps below:

  (a)

  Divide all wells in a play into $i=1,2,\cdots ,n$ well groups, where n is the number of reservoirs. In the Bakken play, for example, there are two main reservoirs. Thus, at this step, we identify two groups of wells,

  Well Sample	Reservoir	Sample Size
  1	Middle Bakken	9860
  2	Three Forks	4818

  and $n=2$ in the Bakken shale.

  (b)

  Further subdivide each of these well groups into $j_i$, $i=1,2,\cdots ,n$, areas with different reservoir qualities. For example, in the Bakken play, the center of the basin has the thickest oil-prolific layer and hence the highest oil production [3,22,28,29,30]. Thus, we have delineated an area at the center of the basin (with maximum oil rate >750 bopd) as core area and the rest as $noncore$ area, see Figure 2 for more details. At this step, we have four well groups that fall into four distinct, static play regions,

  	Reservoir	Region	Sample Size
  1	Middle Bakken	core	5732
  2	Middle Bakken	noncore	4128
  3	Three Forks	core	2672
  4	Three Forks	noncore	2146

  In Bakken, $j_1=2$ and $j_2=2$, that is, there are two reservoir qualities (core and noncore) for each of the reservoirs (Middle Bakken and Three Forks).

  (c)

  Subdivide wells in each of the four regions (two reservoirs × two reservoir qualities each) by time interval classes that encompass significant changes in well completion technology. In the Bakken play, for instance, the newly completed wells have longer lateral lengths, bigger hydraulic fractures and more fracture stages [22,31,32,33,34]. Thus, we classify the wells in each of the four regions by three completion date intervals: [2000–2012], [2013–2016] and [2017–2019]. In the end, we have divided all 14,678 horizontal wells in the Bakken into 12 well groups (4 regions × 3 completion date classes) listed in Table 1. In Bakken, each of the 4 play regions has 3 completion classes, finally yielding 12 static well samples.

• For each well sample, obtain a $P_{50}$ well prototype by fitting a Generalized Extreme Value (GEV) distribution to all qualifying sample wells as follows:

  (a)

  From a given static well sample-k (a region further subdivided by completion dates), consider a dynamic cohort $l_k$ that contains all wells that have at least $l_k$ years on production ($l_k=1,2,\cdots ,t_{{max}_k}$ and $k=1,2,\cdots ,N_{\mathrm{sample}}$). For example: (i) There are 2550 wells in Middle Bakken noncore [2000–2012] group. However, there are only 2540 wells with at least one year on production (the other 10 wells have production records with less than 12 months). Thus, we retain these 2540 wells as cohort-1 of this particular group, see Figure 4. (ii) There are 428 wells in Three Forks core [2000–2012] group. But, only 197 wells have production records of at least seven years. As such, these 197 wells are qualified as cohort-7 of this particular well group, see Figure 5. For detailed GEV fits for all dynamic well cohorts and static well groups in the Bakken play, see Supporting Online Materials-1. For Bakken, $N_{\mathrm{sample}}=12$.

  (b)

  Define a set $X_{l_k}$ that takes values of oil production rate (kbbl/year) from all wells with at least $l_k$ years on production in sample k.

  (c)

  For every $X_{l_k}$, fit a Generalized Extreme Value (GEV) distribution using Equation (A1) and obtain the location parameter, $\mu$, scale parameter, $\sigma$ and shape parameter, $\xi$. (All of these parameters are further subscripted by $l_k$ but we will skip this complication in the notation.)

  (d)

  Calculate ${\overline{x}}_{l_k}$ as the GEV mean of $X_{l_k}$ for each year-$l_k$, $l_k=1,2,\cdots ,t_{{max}_k}$, using Equation (A3).

  (e)

  For each well sample k, obtain the $P_{{50}_k}$ well prototype. We can use the same procedure to obtain other statistic parameters for each well cohort. The $P_{{10}_k}$ and $P_{{90}_k}$ values can be calculated as the ninetieth and tenth percentiles of $X_{l_k}$, respectively. The median and mode can be calculated using Equations (A5) and (A6). by connecting GEV mean values of all years $l_k=1,2,\cdots ,t_{{max}_k}$

  $$P_{{50}_k}=({\overline{x}}_1,{\overline{x}}_2,\cdots ,{\overline{x}}_{t_{{max}_k}}),k=1,2,\cdots ,N_{\mathrm{sample}},\mathrm{kbbl}/\mathrm{year}.$$

  (1)

• Extend $P_{{50}_k}$ (we will subsequently skip the implied subscript k) well prototypes by fitting to them physical scaling curves as follows:

  (a)

  For a each well sample k, calculate the observed cumulative mass produced, $\mathbf{m}$ (ktons) from Equation (A13). Fix $q_{o,\mathrm{ST}}=P_{{50}_k}$ in kbbl/year and $\Delta t_i=1$ year.

  (b)

  Adjust $\tau$ in Equation (A22) or Equation (A7) (with or without exterior flow) and $\mathcal{M}$ in Equation (A11), so that the physical scaling curve matches the observed $\mathbf{m}$. Use $c_t/S_{oi}(p_i-p_f)$ and a values from the average Bakken reservoir properties listed in

  Table A1. Reservoir properties used in scaling oil production in the Middle Bakken and Three Forks (part-1).

  Parameter	Middle Bakken		Three Forks		Data Source
  	SI Units	Field Units	SI Units	Field Units
  Initial pressure, $p_i$	36.8 Mpa	5340 psia	37.1 Mpa	5380 psia	[22]
  Fracture pressure, $p_f$	3.4 Mpa	500 psia	3.4 Mpa	500 psia	[22]
  Connate water saturation, $S_{wc}$	0.57	0.57	0.65	0.65	[22]
  Initial oil saturation, $S_{oi}$	0.43	0.43	0.35	0.35	[22]
  Rock porosity, $\varphi$	0.046	0.046	0.058	0.058	[22]
  Rock permeability, k	4.4 × ${10}^{-17}$ m${}^2$	0.045 md	4.6 × ${10}^{-17}$ m${}^2$	0.047 md	[22]
  Rock compressibility, $c_{\varphi }$	4.3 × ${10}^{-10}$ Pa${}^{-1}$	3.0 × ${10}^{-6}$ psi${}^{-1}$	4.3 × ${10}^{-10}$ Pa${}^{-1}$	3.0 × ${10}^{-6}$ psi${}^{-1}$	[45]
  Water compressibility, $c_w$	4.3 × ${10}^{-10}$ Pa${}^{-1}$	3.0 × ${10}^{-6}$ psi${}^{-1}$	4.3 × ${10}^{-10}$ Pa${}^{-1}$	3.0 × ${10}^{-6}$ psi${}^{-1}$	[45]
  Oil compressibility, $c_o$	1.4 × ${10}^{-9}$ Pa${}^{-1}$	1.0 × ${10}^{-5}$ psi${}^{-1}$	1.4 × ${10}^{-9}$ Pa${}^{-1}$	1.0 × ${10}^{-5}$ psi${}^{-1}$	[45]
  Total compressibility, $c_t$	6.3 × ${10}^{-10}$ Pa${}^{-1}$	9.0 × ${10}^{-6}$ psi${}^{-1}$	5.9 × ${10}^{-10}$ Pa${}^{-1}$	8.5 × ${10}^{-6}$ psi${}^{-1}$
  Oil viscosity, ${\mu }_{o,i}$	3.9 × ${10}^{-5}$ Pa s	0.392 cp	2.8 × ${10}^{-5}$ Pa s	0.276 cp	[46]
  API gravity	42.2${}^{\circ }$API	42.2${}^{\circ }$API	38.7${}^{\circ }$API	38.7${}^{\circ }$API	[46]
  $(c_t/S_{oi})(p_i-p_f)$	0.1014	0.1014	0.1178	0.1178

  and

  Table A2. Reservoir properties used in scaling oil production in the Middle Bakken and Three Forks (part-2).

  Completion	Number of Frac	2d	$\frac{d}{L}$	a
  Date	Stages, N	(ft)		Middle Bakken	Three Forks
  2000–2012	18 [31]	528	0.44	0.050	0.058
  2013–2016	35 [32]	271	0.23	0.026	0.030
  2017–2019	60 [33]	158	0.13	0.015	0.018

  . The matching values of $\tau$ and $\mathcal{M}$ for two both scaling curves and for all well samples are detailed in Table 2.

  (c)

  Get the extended cumulative mass produced, $\widehat{\mathbf{m}}$ (ktons) by multiplying the fitted $\mathcal{M}$ and $\mathrm{RF}(t/\tau )$ from the master curve where t is now calculated each month The benefit of matching $P_{50}$ with a physics-based scaling curve is that we can interpolate and extrapolate production data precisely. Thus, we can change time intervals from years to months and forecast production decades into the future. We recommend to use monthly intervals for precise forecasts that is, $t=\frac{1}{12},\frac{2}{12},\cdots ,50$ years.

  (d)

  Obtain the extended $P_{50}$ well prototypes for well sample k by differencing $\widehat{\mathbf{m}}$ converted to volume

  $${\widehat{P}}_{50}=\frac{\Delta \widehat{\mathbf{m}}}{\Delta t({\rho }_{o,\mathrm{ST}}+R_s{\rho }_{g,\mathrm{ST}})},\mathrm{kbbl}/\mathrm{month}$$

  (2)

• Obtain base forecast of total field oil production for existing wells as follows:

  (a)

  Create a calendar date series with monthly intervals for example, (1/1/2000, 2/1/2000, ⋯, 12/1/2050) and assign $N_{\mathrm{dates}}$ as the length of this series.

  (b)

  Create an empty vector $\mathbf{s}$ with the length of $N_{\mathrm{dates}}$: ($\mathbf{s}=[s_1{s}_2\cdots s_{N_{\mathrm{dates}}}]=[0\cdots 0]$).

  (c)

  For each well, find index $ii$ in the calendar date series that brackets well completion date.

  (d)

  For each well, add to vector $\mathbf{s}$ its corresponding ${\widehat{P}}_{50}$ right-shifted by $ii$.

  $$\mathbf{s}=\mathbf{s}+[\underset{ii\mathrm{zeros}}{\underbrace{0\cdots 0}}{\widehat{P}}_{50}\left(1\right)\cdots {\widehat{P}}_{50}(N_{\mathrm{dates}}-ii+1)]$$

  (3)

  (e)

  Calculate total field oil rate, $\mathbf{q}$ and total field cumulative oil, $\mathbf{Q}$ for existing wells as follows:

  $$\mathbf{q}=\mathbf{s}/(365.25/12)\times {10}^{-3}\mathrm{million}\mathrm{bbl}\mathrm{of}\mathrm{oil}\mathrm{per}\mathrm{day}$$

  (4)

  $$\mathbf{Q}=\left[s_1{s}_1+s_2\cdots \sum _{i=1}^{N_{\mathrm{dates}}}{s}_i\right]\times {10}^{-6}\mathrm{billion}\mathrm{bbl}$$

  (5)

• Calculate infill potential as follows:

  (a)

  Create a one-square-mile fishnet inside the boundary of each area defined in step I(b).

  (b)

  Search all existing wells located inside the boundary.

  (c)

  Calculate wellhead density by counting the number of wells on each of the one-square-mile squares.

  (d)

  Calculate approximate well density following algorithm below.

  • For each well, calculate the number of squares, n, intercepted by the lateral For example, a 5000 ft lateral occupies one square, a 9000 ft lateral occupies two squares and a 14,000 ft lateral occupies three squares because one mile is 5280 ft.
  • Search for the least-occupied n squares in all possible directions.
  • Increase the value of well density by 1 well/mi${}^2$ for every least dense n square found.
  • Repeat points i–iii until all wells in the area of interest are exhausted.

  (e)

  Calculate infill potential by subtracting the calculated well density from the maximum allowable number of wells that still avoid frac hits. For example, in the Bakken play, the tip-to-tip hydraulic fracture length is roughly 1200 ft. 5280 ft/1200 ft ≈ 4. Therefore, the maximum allowable number of wells to avoid frac hits is 4 wells per square mile. Suppose that at some location, the well density already is 3 wells per square mile. Then, the infill potential is $4-3=1$ well per square mile.

• Obtain an infill forecast of total field oil production after adding future wells.

  (a)

  Create monthly drilling schedules for every infill potential and assume a constant drilling rate based on current rig availability, see Figure 9b and SOM-2.

  (b)

  Use the same calendar date series as in step IV(a) with $N_{\mathrm{dates}}$ as the length of date series.

  (c)

  Create an empty vector ${\mathbf{s}}_{\mathrm{f}}$ with the length of $N_{\mathrm{dates}}$ (${\mathbf{s}}_{\mathrm{f}}=[s_{\mathrm{f},1}{s}_{\mathrm{f},2}\cdots s_{\mathrm{f}},\mathrm{N}]=[0\cdots 0]$).

  (d)

  For every drilling schedule, find index $ii$ in the calendar date series that brackets infill schedule date. Let $N_{\mathrm{f}}$ be the number of wells to be drilled on that date.

  (e)

  For every drilling schedule, add to vector ${\mathbf{s}}_{\mathrm{f}}$ its corresponding ${\widehat{P}}_{50}$ after right-shifting it by $ii$ and multiplying by $N_{\mathrm{f}}$.

  $${\mathbf{s}}_{\mathrm{f}}={\mathbf{s}}_{\mathrm{f}}+N_{\mathrm{f}}\times [\underset{ii\mathrm{zeros}}{\underbrace{0\cdots 0}}{\widehat{P}}_{50}\left(1\right)\cdots {\widehat{P}}_{50}(N_{\mathrm{dates}}-ii+1)]$$

  (6)

  (f)

  Calculate total field oil rate, ${\mathbf{q}}_{\mathrm{f}}$ and total field cumulative oil, ${\mathbf{Q}}_{\mathrm{f}}$ after infill drilling as follows

  $${\mathbf{q}}_{\mathrm{f}}=\mathbf{q}+{\mathbf{s}}_{\mathrm{f}}/(365.25/12)\times {10}^{-3}\mathrm{million}\mathrm{bbl}\mathrm{of}\mathrm{oil}\mathrm{per}\mathrm{day}$$

  (7)

  $${\mathbf{Q}}_{\mathrm{f}}=\mathbf{Q}+\left[s_{\mathrm{f},1}{s}_{\mathrm{f},1}+s_{\mathrm{f},2}\cdots \sum _{i=1}^{N_{\mathrm{dates}}}{s}_{\mathrm{f},\mathrm{i}}\right]\times {10}^{-6}\mathrm{in}\mathrm{billion}\mathrm{bbl}$$

  (8)

  where $\mathbf{q}$ and $\mathbf{Q}$ are total field rate and field cumulative from existing wells calculated before with Equations (4) and (5).

5. Conclusions

• We have provided a transparent hybrid method of forecasting oil production at shale basin scale.
• Our statistical approach generates the non-parametric well prototype templates that are used to calibrate our physics-based flow scaling with late-time radial inflow.
• In particular, our average $P_{50}$ well prototypes follow the physics of linear transient flow and are used to calibrate the physics-based scaling extensions to 30 years on production.
• A combination of GEV statistics with physical scaling matches historical production data almost perfectly and gives a smooth, physics-based estimate of future production.
• Our prediction of the Bakken future is optimal in the least square sense [37,38,39,40]; in other words, our prediction is as good as it gets given all data at hand and first order physics of oil recovery in the Bakken shale.
• Regulators may want to consider our approach as a prerequisite to booking reserves in oil shales.
• Newly completed wells have almost the same ultimate recovery as the older ones, despite their much higher initial oil rates.
• Ultimately, we predict that the 14,678 existing wells in the Bakken will produce 5 billion bbl of oil by 2050 (∼340 kbbl/well).
• After drilling additional 4400 new wells at the rate of 120 wells/month, the core area of the Bakken will be drilled out by 2021 and ultimate recovery will be 7 billion barrels of oil.
• With 26,500 more wells to be drilled in the noncore area until 2041, ultimate recovery in the Bakken might be 13 billion barrels of oil but drilling of such scale is unlikely to happen.
• Policy-makers should not of assume that oil boom in the Bakken shale will last decades longer.