Physical Scaling of Oil Production Rates and Ultimate Recovery from All Horizontal Wells in the Bakken Shale

Abstract

A recent study by the Wall Street Journal reveals that the hydrofractured horizontal wells in shales have been producing less than the industrial forecasts with the empirical hyperbolic decline curve analysis (DCA). As an alternative to DCA, we introduce a simple, fast and accurate method of estimating ultimate recovery in oil shales. We adopt a physics-based scaling approach to analyze oil rates and ultimate recovery from 14,888 active horizontal oil wells in the Bakken shale. To predict the Estimated Ultimate Recovery (EUR), we collapse production records from individual horizontal shale oil wells onto two segments of a master curve: (1) We find that cumulative oil production from 4845 wells is still growing linearly with the square root of time; and (2) 6401 wells are already in exponential decline after approximately seven years on production. In addition, 2363 wells have discontinuous production records, because of refracturing or changes in downhole flowing pressure, and are matched with a linear combination of scaling curves superposed in time. The remaining 1279 new wells with less than 12 months on production have too few production records to allow for robust matches. These wells are scaled with the slopes of other comparable wells in the square-root-of-time flow regime. In the end, we predict that total ultimate recovery from all existing horizontal wells in Bakken will be some 4.5 billion barrels of oil. We also find that wells completed in the Middle Bakken formation, in general, produce more oil than those completed in the Upper Three Forks formation. The newly completed longer wells with larger hydrofractures have higher initial production rates, but they decline faster and have EURs similar to the cheaper old wells. There is little correlation among EUR, lateral length, and the number and size of hydrofractures. Therefore, technology may not help much in boosting production of new wells completed in the poor immature areas along the edges of the Williston Basin. Operators and policymakers may use our findings to optimize the possible futures of the Bakken shale and other plays. More importantly, the petroleum industry may adopt our physics-based method as an alternative to the overly optimistic hyperbolic DCA that yields an ‘illusory picture’ of shale oil resources.

1. Introduction

In early 2019, seven major shale plays in the US produced nearly 70% of total US oil production. The Bakken shale is the second largest oil-producing play that in mid-2019 produced 1.45 million bbl of oil per day from 14,888 active horizontal hydrofractured wells. Therefore, it is important to have a reliable predictive model of oil production and ultimate recovery in the Bakken and other shale plays. In shale or tight oil reservoirs, where drainage area and communication between wells are limited [1,2], the estimated ultimate recovery (EUR) per well is the key indicator used to assess fieldwide production potential. The most popular way of predicting EUR from a shale well is the traditional decline curve analysis (DCA) proposed by Arps [3] 75 years ago. This method is favored by the industry, because it only requires production data to forecast future performance of each well and it is easy to implement for many wells. Most published fieldwide predictions of EUR use DCA. Good examples are the Bakken predictions by the US Geological Survey [4,5] and US DOE Energy Information Administration [6].

In his seminal paper, Arps [3] stated that DCA works by extrapolating oil rate observed in a pseudo-steady state (boundary-dominated) flow regime. Unfortunately, because of the vanishing matrix permeability, shale wells do not reach pseudo-state flow regime during their entire existence. Thus, traditional DCA commonly overestimates EURs of shale wells [7,8]. In early 2019, an analysis published by the Wall Street Journal (WSJ) revealed that 2/3 of the projections made by shale companies appeared to be overly optimistic [9,10]. Thousands of shale wells drilled between 2014 and 2017 are pumping on average 10% less oil and gas than their owners’ forecasts. In certain regions, the discrepancy between forecasts and actual production is more than 50%. The total deficit is equivalent to one billion barrels of oil and gas that will not be produced over 30 years or to $30 billion at current oil price.

Several authors have suggested alternatives to the traditional DCA. These are, e.g., the Power Law Exponential Decline [11], the Stretched Exponential Decline [12,13], the Logistic Growth Model [14], the Bayesian probabilistic DCA [15,16], and the Extended Exponential DCA [17]. Empirical models are simple, but they also are mere curve fits, whose parameters have little physical meaning and cannot be generalized. In hyperbolic DCA, for example, the exponent b is adjusted to fit production data, notwithstanding several studies that have tried to elucidate its physical meaning [18,19,20]. Another limitation of most empirical models is their inability to capture the two regimes of oil and gas production in shales that are transient-linear flow at early times and boundary-dominated flow at later times on production [21]. We posit that the physics-based analytic or numeric solutions (such as [22,23]) are more predictive than the empirical decline curves. Unfortunately, many analytical solutions proposed in recent years are unpopular because of their complexity and large quantities of input data.

In this paper, we reintroduce our physics-based scaling approach that is based on an analytical solution of transient oil flow. In several other papers, notably in [14,21,24,25,26,27], we have reviewed the vast literature on this subject and, to avoid repetition, we omit it here. The physical scaling approach was originally developed to predict gas production in the Barnett shale [24,25]. In 2017, we extended our earlier gas production scaling to predict EURs from thousands of oil wells in the Eagle Ford [28]. In this paper, in turn, we extend the Eagle Ford scaling to predict EURs from the 14,888 existing horizontal, hydrofractured black oil wells in the Bakken. We consider the simplest model of oil production consistent with the first order physics and idealized geometry of the extraction process. In principle, the solution of this model depends on many parameters, but in practice—for a given shale play—all but two can be fixed at a field region-average values, resulting in an exact solution of a linearized transient pressure diffusion problem with simplified geometry of the hydrofractured horizontal wells as the boundary condition. At early times, the scaled oil production rate declines as one over the square root of dimensionless time on production, reflecting transient oil flow. Later, oil production rate declines exponentially indicating boundary-dominated flow. In addition, we combine the physical scaling approach with time superposition to resolve erratic production records of the refractured old Bakken wells, and of the wells that were put on artificial lift with time-varying downhole conditions.

To predict EURs, we collapse production records from individual horizontal oil wells onto two different segments of an almost universal scaling curve. The scaling is easy to program and works perfectly for large production data sets. Even though our approach is mechanistic, it preserves the simplicity of traditional DCA by relying only on average reservoir rock and fluid properties and on well-by-well time series of oil production. In a slightly different context, our physical scaling approach has been proven to have the same predictive power as full reservoir simulation, while being several orders of magnitude faster [29]. The EURs calculated from the physical scaling are compared with the latest results published by EIA.

We map water cuts, gas–oil ratios, and calculated EURs of Bakken wells to locate the best producing region in the Williston basin. We also analyze how recent advancements of well completion technologies influence EUR. Operators and policymakers may use our findings to guide sensible development of the Bakken and other shale plays. More importantly, petroleum industry may adopt our simple, fast and accurate method to predict oil production in shales as an alternative to the overly optimistic hyperbolic DCA that yields an ‘illusory picture’ of shale oil resources [9].

2. Materials and Methods

2.1. Physical Scaling

The physical scaling approach extended in this paper was first derived in the paper SPE-187226 [28], and used to predict oil production from thousands of horizontal wells in the Eagle Ford shale play. The scaling solution given in SPE-187226 is based on an assumption that oil is produced from a hydrofractured shale well under a natural depletion mechanism (Section 2.1.1). If a well intervention, most notably refracturing of old wells, occurs, this well’s production cannot be scaled with a single master curve. Therefore, we implement superposition in time to address the erratic production records caused by well refracturing (Section 2.1.2), or putting a well on artificial lift, or changing downhole flowing pressure, or a combination of two or more of these factors.

2.1.1. Physical Scaling of Natural Depletion

Figure 1a illustrates the simple conceptual model used to derive our scaling method. The fracture height, H, is assumed to be the formation thickness; $2L$ is the tip-to-tip hydrofracture length; and $2d$ is the distance between two consecutive hydraulic fractures (or stages). $\mathcal{M}$ is the mass of initial oil in place inside of the stimulated reservoir volume (SRV). In this model, radial flow into the wellbore is neglected, because the hydrofracture permeability is much higher than that of the rock matrix. As a result, we can assume bi-linear flow towards hydraulic fracture planes and write a one-dimensional nonlinear pressure diffusion equation inside SRV. The oil production rate

$$q=-2\frac{k{A}_f}{\mu }{\left(\frac{\partial p}{\partial x}\right)}_{x=-d}$$

(1)

is obtained from the solution of the nonlinear—in general—oil pressure diffusion equation

$$\frac{\partial p}{\partial t}=\alpha \left(p\right)\frac{{\partial }^{2}p}{\partial x^2}$$

(2)

with

$$\begin{array}{cc}\hfill \alpha \left(p\right)& =\frac{k}{\varphi c_t{\mu }_o}\frac{{\rho }_o}{{\rho }_{\mathrm{fluid}}}\hfill \\ \hfill c_t& =S_{oi}{c}_o+S_{wc}{c}_w+c_{\varphi }\hfill \\ \hfill {\rho }_{\mathrm{fluid}}& =\varphi (S_o{\rho }_o+S_w{\rho }_w+S_g{\rho }_g),\frac{\mathrm{kg}}{\mathrm{m}{}^3\mathrm{bulk}}\hfill \end{array}$$

(3)

Here $\alpha$ is the pressure diffusivity coefficient that depends on the reservoir fluid pressure; ${\rho }_{\mathrm{fluid}}$ is the bulk density or reservoir fluid that may contain free gas in the near-hydrofracture regions; and $c_t$ is the total system compressibility near bubble point pressure. For the derivation details, please see [28].

The initial and boundary conditions for this problem are

$$\begin{array}{cc}\hfill p(x,t=0)& =p_i,\mathrm{Initial}\mathrm{pressure}\mathrm{is}\mathrm{uniform}\hfill \\ \hfill p(x\pm 1,t>0)& =p_f,\mathrm{Hydrofracture}\mathrm{pressure}\mathrm{is}\mathrm{constant}\hfill \\ \hfill & \mathrm{during}\mathrm{production}\hfill \end{array}$$

(4)

Briefly, with several assumptions that reflect the simplified physics of oil flowing mostly above the bubble point, the initial and boundary value problem (IBVP) (Equations (2)–(4)) was solved analytically [28], and the following nonlinear least square optimization problem was obtained

$$\begin{array}{cc}\hfill & \underset{\mathcal{M}>0,\tau >0}{min}\mathcal{F}(\mathcal{M},\tau ,\dots )=\hfill \\ \hfill & \sum _{i=1}^{t_{max}}{\left[\left({\rho }_{o,\mathrm{ST}}+R_s{\rho }_{g,\mathrm{ST}}{q}_{o,\mathrm{ST}}\left(t_i\right)\right)\mathsf{\Delta }{t}_i-\frac{c_t\mathcal{M}}{S_{oi}}(p_i-p_f)\left(1-\frac{8}{{\pi }^2}\sum _{n=0}^{\infty }\frac{1}{{(2n+1)}^2}{\mathrm{e}}^{(-{(2n+1)}^2{\pi }^2{t}_i/\left(4\tau \right))}\right)\right]}^2\hfill \end{array}$$

(5)

This problem has two unknowns,

$$\mathcal{M}=2L\left(2dN\right)H{\rho }_{oi}\varphi S_{oi}$$

(6)

and

$$\tau =\frac{d^2}{{\alpha }_i},\mathrm{where}{\alpha }_i=\frac{k}{\varphi c_t{\mu }_{oi}}\frac{{\rho }_{oi}}{{\rho }_{\mathrm{fluid},i}}$$

(7)

where $\mathcal{M}$ is the initial mass of oil in place inside the stimulated reservoir volume; $2dN$ is the total length of horizontal well; and $\tau$ is the characteristic time of pressure interference, or the square of interfracture half-distance, d, divided by the hydraulic diffusivity, ${\alpha }_i$, at initial reservoir conditions. The parameter $\tau$ is the characteristic time needed for the depletion of oil pressure to diffuse to the midplane between two consecutive hydrofractures.

The complex, highly nonlinear optimization problem in Equation (5) is difficult to solve numerically. Thus, a simpler but equivalent “physical scaling curve method” was introduced and illustrated in Figure 1 (right). This method consists of two steps. First, a master curve is constructed using Equation (8) with dimensionless time, $\tilde{t}=\frac{t}{\tau }$, as the x-axis and dimensionless recovery factor, RF, as the y-axis. Second, for each well, its recorded cumulative produced mass is matched to the master curve by scaling the x-axis with $\tau$ and y-axis with $\mathcal{M}$. Please note that these two steps are equivalent to solving the nonlinear optimization problem in Equation (5). The left hand side of Equation (5), ${\rho }_{o,\mathrm{ST}}+R_s{\rho }_{g,\mathrm{ST}}{q}_{o,\mathrm{ST}}\left(t_i\right)$, is the monthly mass production, $N_p$; and the right hand side is the recovery factor in Equation (8) multiplied by the accessible mass in place, $\mathcal{M}$ (remember that $\tilde{t}=\frac{t}{\tau }$).

$$\mathrm{RF}\left(\tilde{t}\right)=C\left(1-\frac{8}{{\pi }^2}\sum _{n=0}^{\infty }\frac{1}{{(2n+1)}^2}{\mathrm{e}}^{(-{(2n+1)}^2{\pi }^2\tilde{t}/4)}\right)$$

(8)

with

$$C=\frac{c_t}{S_{oi}}(p_i-p_f)$$

(9)

being the ultimate oil recovery.

To make our scaling method even more ‘engineer-friendly’, we provide the best closed-form approximation of the exact solution in Equation (8):

$$\mathrm{RF}\left(\tilde{t}\right)=C(1-exp[1-exp\left(\lambda {\tilde{t}}^{\beta }\right)])$$

(10)

For the Bakken, the optimal values of constants $\lambda$ and $\beta$ are 1.3 and 0.55, respectively. The universal mathematical constant $\lambda \sim 1.3$ was obtained by Conway in his lost Cosmological Theorem, which asserts that after a certain number of steps from the Big Bang—the beginning of a Conway sequence of integers [30,31]—all the exotic elements, all the things that are not compounds of common atoms, of digit 1 (hydrogen) for example, disappear and only common atoms are left, also see [32] pp. 452–455. The power exponent, $\beta$, was obtained from a fit. Equation (10) is much easier to calculate than the summation of the infinite series in Equation (8), and may be computed with Excel${}^{\circledR }$.

Now, the linearized oil pressure diffusion equation can be solved painlessly, well-by-well, and the resulting set of fit parameters, $\{\tau ,\mathcal{M}\}$, provides key insights into the physics of production of the mostly saturated shale oil from thousands of wells. Finally, each EUR is calculated as the maximum recovery factor value from the master curve, $\mathrm{RF}\left({\tilde{t}}_{max}\right)$, multiplied by the mass of oil in place, $\mathcal{M}$.

$$\begin{array}{cc}\hfill \mathrm{EUR}& =\mathrm{RF}\left({\tilde{t}}_{max}\right)\mathcal{M}=C\mathcal{M}\hfill \end{array}$$

(11)

2.1.2. Physical Scaling of Well Refracturing

In many cases, the first hydraulic fracturing job is unsuccessful and results in a poor quality well that is a candidate for refracturing at a later time. Refracturing contributes to sudden jumps in production rate. The refracturing process is illustrated as an increase of the SRV in Figure 2a. As before, $\mathcal{M}$ is the initial mass of oil in place in the SRV. Before refracturing, the SRV box contains ${\mathcal{M}}_0$ ktons of oil in place. At time $t_1$, the first refracturing job is completed, and the mass in place jumps to ${\mathcal{M}}_1$. At $t_2$, the second refracturing job is performed, and an even bigger ${\mathcal{M}}_2$ is accessed. The refracturing process increases the mass $\mathcal{M}$ used to rescale cumulative oil production. By superposition theorem, this oil production is a summation of the production before refracturing that accessed the original mass in place, ${\mathcal{M}}_0$, and secondary production from imaginary wells at times $t_j$, $j=1,2,\dots$ that access the incremental masses in place, (${\mathcal{M}}_j-{\mathcal{M}}_{j-1}$).

$$\begin{array}{cc}\hfill & \underset{\mathcal{M}>0,\tau >0,a_j\ge 0}{min}\mathcal{F}(\mathcal{M},\tau ,a_j,\dots )=\sum _{i=1}^{t_{max}}[\left({\rho }_{o,\mathrm{ST}}+R_s{\rho }_{g,\mathrm{ST}}{q}_{o,\mathrm{ST}}\left(t_i\right)\right)\mathsf{\Delta }{t}_i-\hfill \\ \hfill & \sum _{j=0}^{n_{\mathrm{jump}}}\left(a_j{\mathbb{H}}_j\frac{c_t{\mathcal{M}}_0}{S_{oi}}(p_i-p_f)\left(1-\frac{8}{{\pi }^2}\sum _{n=0}^{\infty }\frac{1}{{(2n+1)}^2}{\mathrm{e}}^{(-{(2n+1)}^2{\pi }^2{t}_i-t_j/\left(4\tau \right))}\right)\right){]}^2\hfill \end{array}$$

(12)

The production jumps from refracturing cannot be scaled with the standard master curve shown in Figure 1. Thus, we combine the physical scaling solution with superposition and obtain Equation (12) that remedies this difficulty. If primary production before refracturing has the recovery constant C in Equation (13), then each refracturing event j at time $t_j$ activates an imaginary well with the recovery constant $C{a}_j$. Here $a_j$’s are the new fitting parameters that scale effectiveness of each refracturing job relative to the initial mass in place, Equation (14). ${\mathbb{H}}_j$ is the Heaviside step function that handles time delays in the consecutive refracturing jobs.

$$C=\frac{c_t}{S_{oi}}(p_i-p_f)$$

(13)

$$a_j=\frac{{\mathcal{M}}_j-{\mathcal{M}}_{j-1}}{{\mathcal{M}}_0}$$

(14)

$${\mathbb{H}}_j=\left\{\begin{array}{cc}0,\hfill & t_i-t_j<0,\hfill \\ 1,\hfill & t_i-t_j\ge 0,\hfill \end{array}\right.$$

(15)

The total master curve ${\mathrm{RF}}_T$ is the summation of all individual master curves, ${\mathrm{RF}}_j$, see Equations (16) and (17). An example of matching the scaled discontinuous historical production is shown in Figure 2b.

$${\mathrm{RF}}_T\left(\tilde{t}\right)=\sum _{j=0}^{n_{\mathrm{jump}}}{\mathrm{RF}}_j\left(\tilde{t}\right)$$

(16)

$$\begin{array}{cc}\hfill & {\mathrm{RF}}_j\left(\tilde{t}\right)=a_j{\mathbb{H}}_{j}C\left(1-\frac{8}{{\pi }^2}\sum _{n=0}^{\infty }\frac{1}{{(2n+1)}^2}{\mathrm{e}}^{(-{(2n+1)}^2{\pi }^2\tilde{t}/4)}\right),\mathrm{where}\tilde{t}=\frac{t_i-t_j}{\tau }\hfill \end{array}$$

(17)

Finally, one calculates EUR by multiplying the maximum total recovery factor from the modified master curve, ${\mathrm{RF}}_T\left({\tilde{t}}_{max}\right)$, with the initial mass of oil in place before fracturing, ${\mathcal{M}}_0$, Equation (18).

$$\begin{array}{cc}\hfill {\mathrm{EUR}}_T& ={\mathrm{RF}}_T\left({\tilde{t}}_{max}\right){\mathcal{M}}_0=C\left(1+\sum _{j=1}^{n_{\mathrm{jump}}}{a}_j\right){\mathcal{M}}_0\hfill \end{array}$$

(18)

2.2. Data Collection and Scaling Procedure

We have mined public well records for the Bakken shale from the DrillingInfo (now Enverus) 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 with a premium subscription. In DrillingInfo, we have selected only 14,888 horizontal oil wells that were completed after 1 January 2000, in the Middle Bakken and Three Forks formations. The raw production data from DrillingInfo were imported into an integrated MATLAB${}^{\circledR }$ software package we designed to perform data cleanup, consolidation, and unit conversions. The clean production data for each well consist of a vector of elapsed time on production (consecutive months with nonzero oil volume records) and a vector of cumulative oil production in kilotons.

Our scaling procedure of obtaining EUR for each well is as follows:

• Exclude all newly completed wells with less than 12 months of production.
• For each remaining well, plot its cumulative production vs. square root of time on production. Classify these wells as:

  (a)

  Non-interfering wells, if the plot shows a straight line (with, e.g., ${\mathrm{R}}^2\ge 0.9$).

  (b)

  Interfering wells, if the plot shows a deviation from a straight line (with, e.g., ${\mathrm{R}}^2<0.9$)

  (c)

  Refrac wells, if production jumps are detected.

• For each non-interfering well, scale its cumulative production by $\mathcal{K}\sqrt{t_{max}}$, on the $y-$axis and the elapsed time on production by $t_{max}$ on the $x-$axis to match the line $f\left(\tilde{t}\right)=\sqrt{\tilde{t}}$. To predict EUR, assume that deviation from the line starts at $\tau =2t_{max}$. Thus, the corresponding $\mathcal{M}$ can be calculated as $\mathcal{M}=\left(\mathcal{K}\sqrt{t_{max}}\right)/\mathrm{RF}\left(\frac{1}{2}\right)$, where $\mathrm{RF}\left(\frac{1}{2}\right)$ is calculated from Equation (8). Finally, $\mathrm{EUR}$ is calculated from Equation (11).
• For each interfering well, scale its cumulative production records to the master curve in Equation (8) using two scaling parameters, $\tau$ and $\mathcal{M}$. The corresponding EURs are calculated from Equation (11).
• For each refrac well, scale its cumulative production records to the master curve in Equation (16), and obtain optimum values of the scaling parameters $\tau$, ${\mathcal{M}}_0$, and $a_j$, $j=1,2,\dots$. The corresponding EURs are calculated from Equation (18).
• For each newly completed well, calculate its EUR using expected values of $\tau$ and $\mathcal{M}$ from comparable interfering wells that were completed between 2017 and 2019. Use Equation (11).

Notice that the EURs calculated from Equations (11) and (18) are in mass units. We use oil density to convert these EURs into thousands or millions of barrels (kbbl or Mbbl).

3. Results and Discussion

A quick scoping analysis shows the distribution of oil rates from all Bakken wells after—in this example—one year on production, see Figure 3. This distribution is positively skewed, meaning that most Bakken wells have low oil rates and the best wells with the significantly higher oil rates in the long tail are rare. The Generalized Extreme Value (GEV) distribution is most suitable to model the empirical oil rate distribution in the Bakken and other shales [27,33,34].

Using the GEV statistics of well populations, we calculate the high ($P_{10}$), mean ($P_{50}$), and low ($P_{90}$), oil rate values from the cumulative distribution function (CDF). The $P_{50}$ well is the expected value of the distribution; the $P_{10}$ well is exceeded by only 10% of wells; and the $P_{90}$ well is exceeded by 90% of wells. Figure 4 shows the daily oil rates and cumulative oil production for all existing horizontal wells in the Middle Bakken and the Upper Three Forks. Notice that the Three Forks wells have fewer years on production because they were drilled after 2010. The $P_{10}$, $P_{50}$, and $P_{90}$ wells were calculated with the method discussed elsewhere [27,34] and in the caption of Figure 3. The expected ($P_{50}$) maximum oil rate and cumulative oil are about 500 bbl/d and 200,000 bbl of oil for both Bakken and Three Forks.

In addition, many production jumps are observed. Each jump is a sudden discontinuity of oil rate that is sometimes larger than the initial oil rate in the same well, and an abrupt brake in cumulative production curve. These discontinuities may be caused by refracturing of old wells [35,36,37,38,39] or by changes of downhole well flowing pressure [40,41,42,43].

3.1. Physical Scaling Matches

We have scaled oil production from all 14,888 horizontal oil wells in the Bakken following the procedure detailed in Section 2.2. The reservoir properties used in scaling oil production are listed in

Table A2. Reservoir properties used in scaling oil production in the Middle Bakken and Upper Three Forks.

Parameter	Middle Bakken		Upper Three Forks		Data Source
	SI Units	Field Units	SI Units	Field Units
Horizontal well length, $L_w$	2900 m	9500 ft	2900 m	9500 ft	DrillingInfo
Number of fracture stages, N	30	30	30	30	DrillingInfo
Fracture height, H	10 m	33 ft	12 m	40 ft	well log
Tip-to-tip fracture length, $2L$	360 m	1200 ft	360 m	1200 ft	DrillingInfo
Reservoir temperature, T	113 $°$C	237 $°$F	115 $°$C	239 $°$F	[65]
Initial pressure, $p_i$	36.8 Mpa	5340 psia	37.1 Mpa	5380 psia	well log
Saturation pressure, $p_b$	17.4 Mpa	2530 psia	12.1 Mpa	1753 psia	[65]
Fracture pressure, $p_f$	3.4 Mpa	500 psia	3.4 Mpa	500 psia	DrillingInfo
Connate water saturation, $S_{wc}$	0.57	0.57	0.65	0.65	well log
Initial oil saturation, $S_{oi}$	0.43	0.43	0.35	0.35	well log
Rock porosity, $\varphi$	0.046	0.046	0.058	0.058	well log
Rock permeability, k	4.4 $\times {10}^{-17}$ m${}^2$	0.045 md	4.6 $\times {10}^{-17}$ m${}^2$	0.047 md	well log
Rock compressibility, $c_{\varphi }$	4.3 $\times {10}^{-10}$ Pa${}^{-1}$	3.0 microsip	4.3 $\times {10}^{-10}$ Pa${}^{-1}$	3 microsip	[66]
Water compressibility, $c_w$	4.3 $\times {10}^{-10}$ Pa${}^{-1}$	3 microsip	4.3 $\times {10}^{-10}$ Pa${}^{-1}$	3 microsip	[66]
Oil compressibility, $c_o$	1.4 $\times {10}^{-9}$ Pa${}^{-1}$	1 microsip	1.4 $\times {10}^{-9}$ Pa${}^{-1}$	1 microsip	[66]
Oil viscosity, ${\mu }_{oi}$	3.9 $\times {10}^{-5}$ Pa s	0.392 cp	2.8 $\times {10}^{-5}$ Pa s	0.276 cp	[65]
Oil formation volume factor, $B_{oi}$	1.61 m${}^3$/sm${}^3$	1.61 rbbl/stb	1.48 m${}^3$/sm${}^3$	1.48 rbbl/stb	[65]
API gravity	42 $°$ API	42 $°$ API	39 $°$ API	39 $°$ API	[65]
GOR	1.48 sm${}^3$/sm${}^3$	125 scf/stb	110 sm${}^3$/sm${}^3$	620 scf/stb	[65]

. The key scaling results are shown in Figure 5, Figure 6 and Figure 7.

At the time of this analysis, 4845 Bakken wells, 3349 in the Middle Bakken and 1496 in the Upper Three Forks, were classified as non-interfering wells, because their cumulative production was still growing linearly versus the square root of time. The hyperbolic decline ($\propto 1/\sqrt{\tilde{t}}$) of oil production rate from these wells is caused by the purely transient diffusion of oil pressure in the reservoir. Notice that the x-axis in Figure 5 is scaled with the maximum time on production, $t_{max}$ instead of $\tau$ and the y-axis is scaled with $\mathcal{K}\sqrt{t_{max}}$ to match the $f\left(\tilde{t}\right)=\sqrt{\tilde{t}}$ master curve that intersects point (1,1). In the Middle Bakken, the expected values of $t_{max}$ and $\mathcal{K}\sqrt{t_{max}}$ are 75 months and 30 ktons, respectively, while in Three Forks these values are 65 months and 27 ktons. For the non-interfering wells, we assume arbitrarily that they start to decline exponentially (go under the square-root-of-time linear trend) at $\tau =2\times t_{max}$.

The interfering wells exhibit clear pressure interference between consecutive hydrofractures and are in pseudo-steady state flow regime. For each well k, this interference already occurred at the elapsed time on production $t={\tau }_k$ and, at the time of this analysis, each oil production rate was declining exponentially. For these wells, the recovery factors start bending down at the dimensionless time equal $0.5<0.64$ in [24] because oil is less compressible than gas. The interfering wells scale with the master curve in Figure 1. Of the 6401 interfering wells in the Bakken, 4245 were in the Middle Bakken and 2156 in the Upper Three Forks (Figure 6). The expected values of $\tau$ and $\mathcal{M}$ were 100 months and 420 ktons for the Middle Bakken, and 90 months and 320 ktons for the Upper Three Forks.

Many wells in the Bakken (2363) cannot be scaled to a non-interfering or interfering segment of the scaling curve. The poor matches are caused by multiple production jumps that may be due to refracturing of old wells [35,36,37,38,39]. We note that these jumps could also result from artificial lift and step decreases of downhole well pressure [40,41,42,43]. Each of these wells is scaled individually with a modification of the full physical scaling curve coupled with superposition (illustrated in Figure 2). Four examples of “refrac wells” are shown in Figure 7. The wells in panels a and c have been completed in the Middle Bakken formation, while the those in panels b and d in the Upper Three Forks formation. As we can see, the modified scaling curve matches the discontinuous production records. The magnitude of each a in Equation (14) may reflect effectiveness of a specific refracturing job. For all “refrac wells,” the expected value of a is about 0.6, a little better than 50%.

The remaining 1279 newly completed wells had less than 12 months on production at the time of this analysis. The oil production rates from these wells are usually unstable due to backflow of fracturing fluids, excessive choking and/or lack of takeaway capacity [24,25,44,45,46,47,48].

During the first three months on production in each well, the produced cumulative oil volume is commonly below the square-root-of-time line due to the unloading of SRV from fracturing fluids, whose flow blocks oil production, [24,25]. If a well does not have enough of stable production history, its scaling will underestimate both $\tau$ and $\mathcal{M}$. Therefore, we have decided to use the average $\tau$ and $\mathcal{M}$ values from scaling the recently completed (in 2017–2018) interfering wells with $\tau =45$ months and $\mathcal{M}=495$ ktons for the Middle Bakken and 43 months and 390 ktons for the Upper Three Forks. Notice that the magnitude of $\tau$ is lower and that of $\mathcal{M}$ is higher for the new wells. These shifts are caused by the improvements of well completion technology over the last decade [49,50,51,52,53,54,55], i.e., by bigger hydrofractures, fewer clusters per stage, more fracture stages, etc.

Table A3. Summary of pressure interference time, $\tau$, and mass of oil in place, $\mathcal{M}$, for four different well classes in the Middle Bakken and Upper Three Forks.

Well Class	Middle Bakken							Upper Three Forks
	$\mathit{\tau }$ (months)			$\mathcal{M}$ (ktons)			Number	$\mathit{\tau }$ (months)			$\mathcal{M}$ (ktons)			Number
	${\mathit{P}}_{10}$	${\mathit{P}}_{50}$	${\mathit{P}}_{90}$	${\mathit{P}}_{10}$	${\mathit{P}}_{50}$	${\mathit{P}}_{90}$	of Wells	${\mathit{P}}_{10}$	${\mathit{P}}_{50}$	${\mathit{P}}_{90}$	${\mathit{P}}_{10}$	${\mathit{P}}_{50}$	${\mathit{P}}_{90}$	of Wells
Interfering	170	100	50	720	420	190	4245	140	90	40	540	320	150	2156
Non-Interfering	250	220	180	860	480	250	3349	250	210	160	680	400	170	1496
Refracs	230	150	70	1210	660	230	1549	230	140	70	980	540	200	814
Newly completed	60	50	30	800	520	270	751	60	40	30	620	400	200	528
All wells	200	150	90	850	490	220	9894	180	130	80	660	390	170	4994

details the interference times, $\tau$, and the masses of oil in place, $\mathcal{M}$, for the four well classes: (1) interfering, (2) non-interfering, (3) refracs, and (4) newly completed wells in the Middle Bakken and Upper Three Forks. The surface location maps of these four well classes are shown in Figure 8 for both reservoirs. For the Middle Bakken, wells in the first three classes spread extensively from the center of the Williston Basin in North Dakota to the Elm Coulee Field in Montana, while in the Three Forks few wells were drilled in Montana. Interestingly, the newly completed wells are located mostly in the central part of the Williston Basin in North Dakota. There were some new wells completed in Montana; however, they became inactive within the first few months on production. In general, the magnitudes of $\tau$ and $\mathcal{M}$ in the Upper Three Forks are lower than those in the Middle Bakken, indicating that the Three Forks wells decline faster and have lower EURs than those in the Middle Bakken. The refrac wells have the highest values of $\mathcal{M}$, because of the incremental stimulated reservoir volumes due to refracs. The non-interfering wells have slightly higher EURs than the interfering wells, because of the larger $\tau$ and $\mathcal{M}$ values. As mentioned before, the newly completed wells have shorter $\tau$s and larger $\mathcal{M}$s than all other well classes. Each of these wells starts from a very high initial production rate that declines very fast.

3.2. EUR Predictions

To predict the fieldwide EUR, all individual scaling curves for the 14,888 existing wells in the Bakken are summed up vs. calendar time (Figure 9). For each well class, we note the high fidelity of the physical scaling approach in retracing old production and—by implication—forecasting future production of the existing wells. From Figure 9 (left), we infer that the second production peak in the last three years is not only the contribution from the newly completed wells that are highly productive, but also from the old wells that have been refractured or pumped off recently. From Figure 9 (right), it follows that the expected estimated ultimate recovery from all 14,888 wells in the Bakken play is about 715 million m${}^3$ (4.5 billion bbl) of oil. The EURs for the interfering and non-interfering well classes are almost the same at about 240 million m${}^3$ (1.5 billion bbl) each, as the respective well numbers are close: 6401 and 4845. The refrac wells have a larger average EUR per well with total EUR of 160 million m${}^3$ (1 billion bbl) from only 2363 wells. The recent 1279 wells also have a higher EUR per well with the total EUR of 80 million m${}^3$ (0.5 billion bbl). One can argue that success of the newly completed wells in yielding the highest EURs is not only due to the significant improvements in well completion techniques, but also to the favorable drilling location selections across the most prolific area in the Bakken, near the center of the Williston Basin (recall Figure 8).

Figure 10 shows that our physical scaling method and the hyperbolic decline curve analysis (DCA) give different estimates of oil production from the Bakken shale. To solve for the matching parameters of hyperbolic DCA, initial oil rate, $q_i$; initial decline rate, $D_i$; and the hyperbolic exponent, b, we used the algorithm described in the Reservoir Engineering Handbook (Fourth Edition) by Tarek Ahmed [56]. In addition, we wrote a MATLAB^® program to automate the hyperbolic DCA matches of 14,888 Bakken oil wells. Figure 10a,b are examples of the respective matches of oil production from well -25-083-XXXXX. The rate is plotted on a semilog scale and there seem to be no large discrepancies between the physical scaling (blue curve) and hyperbolic DCA (red curve) matches of the past oil rate (black line). However, when we integrate the rate estimates from both methods to obtain cumulative production, the hyperbolic DCA fails to trace past production and gives an overly optimistic cumulative oil of 64,000 m${}^3$ (400 kbbl) at the last production time recorded. In contrast, the physical scaling yields 55,600 m${}^3$ (350 kbbl) matching exactly historical data. We also quality-checked both methods for 11,246 non-refractured wells in the Bakken, see Figure 10c,d. For each well, each point in these cross-plots quantifies agreement between the actual annual rate (the volume of oil produced during 12 months on production in kbbl/year), in year $1,2,\dots$, and the corresponding values predicted with both methods. There are 61,230 points in each of the two cross-plots. Ideally, these points would follow the diagonal line, $y=x$. For DCA, best data fit, $y=1.12x+1.14$, (the blue line) is biased towards excessive predicted production volumes. The linear model explains the data almost perfectly and its 95% confidence interval (CI) is not plotted. The two dashed blue lines denote the 95% CI band that a new data point in the same set of wells would fall inside of this band. The physical scaling model is not biased, $y=0.99x-0.13$, and its 95% CI is narrower. Finally, we stacked up all rate estimates from the 14,888 oil wells in the Bakken to obtain total field rate shown in Figure 10e. At the field scale, the systematic upward bias of hyperbolic DCA is amplified, thus, the red line fails to track the black line. In contrast, the stack of individual physical scaling matches (blue line) traces historical production. The total field cumulative is shown in Figure 10f. In the end, there is a big discrepancy between both methods of predicting total field EUR from the same field production data set. The physical scaling predicts the EUR from the 14,888 existing wells at 715 million m${}^3$ (4.5 billion bbl), while the hyperbolic DCA is overly optimistic, giving a 45% higher estimate of 1 billion m${}^3$ (6.5 billion bbl).

The current model is a mechanistic scaling of well production in the Bakken shale based on the first order physics of production mechanisms in hydrofractured horizontal wells and on geology. It served as the physical foundation of the statistical prediction [34] of ultimate production in the Bakken based on the Generalized Extreme Value (GEV) statistics. The current scaling model provides an unbiased match of all historical data from all wells in the Bakken. The width of the 95% CI in Figure 10d reflects the natural variability of reservoir and well properties that influence the model predictions around their expected values listed in Table A2. The 61,230 points that make the comparison in Figure 10d demonstrate that the stochastic variability of reservoir and well properties is rather small and the model does an excellent job finding the unbiased values of original oil in place $\mathcal{M}$ and pressure interference time $\tau$ for each of the 14,888 wells in the Bakken.

To compare our results with EUR predictions from the US DOE Energy Information Administration [6], we have remapped each Bakken assessment unit onto counties and calculated the expected value of EUR for each county. The details are tabulated in

Table A4. Bakken EUR Summary.

AU	State	County	EIA 2019				Physical Scaling
			Bakken		Three Forks		Bakken				Three Forks
			EUR	Potential	EUR	Potential	P10	P50	P90	Existing	P10	P50	P90	Existing
			(Mb/well)	wells	(Mb/well)	wells	(Mb/well)	(Mb/well)	(Mb/well)	wells	(Mb/well)	(Mb/well)	(Mb/well)	wells
Central Basin	MT	Daniels	60	189	177	35,870	-	-	-	-	-	-	-	-
Central Basin	MT	McCone	60	528			196	196	196	1	-	-	-	-
Central Basin	MT	Richland	232	1083			318	210	103	102	234	118	44	3
Central Basin	MT	Roosevelt	263	4985			330	212	103	138	228	129	38	11
Central Basin	MT	Sheridan	49	753			359	163	37	4	64	59	53	2
Central Basin	ND	Divide	241	12			-	-	-	-	99	-307	27	3
Central Basin	ND	Dunn	184	73			459	301	156	41	375	254	160	14
Central Basin	ND	McKenzie	263	3749			465	298	158	1540	399	250	125	669
Central Basin	ND	Williams	262	3070			463	300	162	1261	376	250	143	524
Eastern Transitional	ND	Burke	7	2706			267	153	62	117	235	148	81	24
Eastern Transitional	ND	Divide	140	658			237	140	64	96	311	182	78	206
Eastern Transitional	ND	Dunn	423	1099			620	354	125	248	554	335	208	111
Eastern Transitional	ND	Hettinger	169	7			-	-	-	-	-	-	-	-
Eastern Transitional	ND	McLean	623	245			523	310	144	40	1198	501	193	4
Eastern Transitional	ND	Mercer	13	144			-	-	-	-	-	-	-	-
Eastern Transitional	ND	Mountrail	232	2679			683	388	155	1163	464	274	114	201
Eastern Transitional	ND	Stark	169	371			177	177	177	1	125	21	11	2
Eastern Transitional	ND	Ward	80	111			209	209	209	1	-	-	-	-
Elm Coulee–Billings Nose	MT	McCone	80	116			-	-	-	-	-	-	-	-
Elm Coulee–Billings Nose	MT	Richland	183	3421			423	230	75	935	112	64	15	2
Elm Coulee–Billings Nose	ND	Billings	60	828			263	192	17	10	265	149	59	37
Elm Coulee–Billings Nose	ND	Golden Valley	476	130			70	70	70	1	502	337	172	12
Elm Coulee–Billings Nose	ND	McKenzie	184	2449			269	157	62	79	291	134	9	8
Nesson–Little Knife	ND	Billings	167	586			268	135	31	49	280	165	76	82
Nesson–Little Knife	ND	Burke	188	680			258	174	95	65	265	173	91	68
Nesson–Little Knife	ND	Divide	115	603			246	159	82	128	236	151	77	149
Nesson–Little Knife	ND	Dunn	324	2685			615	358	185	1258	569	345	181	608
Nesson–Little Knife	ND	Hettinger	223	106			-	-	-	-	-	-	-	-
Nesson–Little Knife	ND	McKenzie	329	1520			651	397	200	1072	610	362	164	1006
Nesson–Little Knife	ND	Mountrail	312	530			545	326	147	975	503	304	146	655
Nesson–Little Knife	ND	Slope	120	167			-	-	-	-	-	-	-	-
Nesson–Little Knife	ND	Stark	156	2164			199	167	134	2	343	203	93	219
Nesson–Little Knife	ND	Williams	178	1828			365	213	96	410	349	211	100	243
Northwest–Transitional	MT	Daniels	82	2584			-	-	-	-	-	-	-	-
Northwest–Transitional	MT	McCone	82	161			-	-	-	-	-	-	-	-
Northwest–Transitional	MT	Roosevelt	82	1312			-	-	-	-	-	-	-	-
Northwest–Transitional	MT	Sheridan	50	2857			195	100	24	27	67	67	67	1
Northwest–Transitional	MT	Valley	1	1005			-	-	-	-	-	-	-	-
Northwest–Transitional	ND	Divide	180	601			248	158	87	41	239	153	77	110
Northwest–Transitional	ND	Williams	212	669			311	185	81	89	317	224	85	20
AVERAGE EUR (thousand bbl)			189		177		515	309	146		456	278	136
TOTAL WELL				49,464		35,870				9894				4994
TOTAL EUR (billion bbl)			9.3		6.3		5.1	3.1	1.4		2.3	1.4	0.7

. We predict the EUR from all 14,888 active horizontal oil wells in the Bakken to be 715 million m${}^3$ (4.5 billion barrels of oil) or ∼48 thousand m${}^3$/well (300 kbbl/well), with 493 million m${}^3$ (3.1 billion barrels) coming from 9894 wells in the Middle Bakken and 222 million m${}^3$ (1.4 billion bbl) coming from 4994 wells in the Upper Three Forks. The results in Table A4 are also displayed as stacked bar graphs for the Middle Bakken and the Upper Three Forks (Figure 11). The stacked bars show the numbers of existing wells from the DrillingInfo dataset in 2019, and the potential wells from [6] in each Bakken Total Petroleum System (TPS) assessment unit mapped onto the corresponding counties in the Middle Bakken and Upper Three Forks Formations. The lines show the corresponding EUR values from the physical scaling and the EIA statistical predictions. The $P_{10}$, $P_{50}$, and $P_{90}$ wells were obtained by applying the generalized extreme value (GEV) distribution [27,34] for each set of EURs from the scaling curve method in each sub-region as demonstrated before in Figure 3.

In most cases, where there are enough wells per sample, the expected $P_{50}$ EURs from physical scaling agree with the EUR estimates by EIA. One can observe that most of the best regions with the high EUR values have been extensively drilled leaving no or few potential wells to infill. Mountrail-NL is an extreme example of over-drilling in the best region of the Bakken, where the number of existing wells exceeds the number of EIA potential wells by 40% (the green bar is taller than the gray bar). In another paper [34], we have calculated well density of all regions in the Bakken, revealing a surprising result that the core area of the Bakken will be fully drilled by 2021. In contrast, most of the 85,334 potential wells proposed by EIA will be in the poorest regions of the Bakken. Figure 12 maps the EUR results from physical scaling approach for 14,888 horizontal oil wells in the Bakken. Notice that this map resembles more the water cut map (Figure A3), rather than the GOR map (Figure A3).

In Figure 13a, we show that EUR is negatively correlated with water cut ($1-\mathrm{oil}\mathrm{cut}$). On the other hand, there is little correlation between EUR and lateral length, as depicted in Figure 13b. Doubling lateral length from 1.5 to 3 km (5000 to 10,000 ft) does not double average EUR that increases from 38 to 49 thousands of m${}^3$ (240 to 310 kbbl). This finding may also suggest that better technology cannot increase EUR if a well is drilled in a low-productivity region of the Bakken.

Figure 14 shows the average values of pressure interference time, $\tau$, lateral length, gas–oil ratio, and water cut for each six-month time interval in calendar dates. The value of $\tau$ decreases significantly over time meaning that newer wells pressure-interfere and decline faster than the old ones. Equation (7) shows that $\tau$ is directly proportional to the square of half-distance, d, between two consecutive hydrofractures. This means that halving the interfracture distance decreases $\tau$ four-fold. Historically, the number of hydrofracture stages in the Bakken has increased over time from 8 stages in 2007 [52], to 18 stages in 2009 [53], 35 stages in 2016 [54], and to as many as 60 stages in 2019 [55]. The lateral length has also doubled from 1.5 km (5000 ft) in 2005 to 3 km (10,000 ft) in 2015.

In summary, the new completion techniques in the Bakken shale rely on more of bigger hydraulic fractures and on longer laterals [49,50,51]. These new completions have significant side-effects. Cumulative gas–oil ratio has been increasing in both the Middle Bakken and Three Forks. The likely culprits are reservoir depressuring by hydrofractures that are too close and the increasingly aggressive artificial lift. Degassing the low permeability and porosity reservoirs in the Bakken spells trouble for future production, because solution- and free gas are the main driving mechanisms of primary depletion. Water cut has also increased with time revealing another negative side-effect of the recent Bakken field development. The increasing water cut might be caused by drilling new low-productivity wells in the immature regions of Bakken [57,58,59], and/or by the longer and taller hydraulic fractures that break into adjacent water-bearing formations, the Lodgepole and/or Lower Three Forks [60,61,62].

4. Conclusions

We have delivered a comprehensive physics-based scaling of all horizontal hydrofractured wells completed in the Middle Bakken and Upper Three Forks between 1 January 2000, and 1 May 2019. In particular:

• The current 14,888 active oil wells in the Bakken shale will ultimately produce 715 million m${}^3$ (4.5 billion barrels) of oil, with 493 million m${}^3$ (3.1 billion barrels) from 9894 wells in the Middle Bakken and 222 million m${}^3$ (1.4 billion barrels) from 4994 wells in the Upper Three Forks.
• In general, wells completed in the Middle Bakken produce more oil than those in the Upper Three Forks due to: (1) lower water saturation and water cut, (2) slower decline rate (longer pressure interference times, $\tau$), and (3) higher initial oil in place (larger $\mathcal{M}$).
• Newly completed wells start from very high initial oil rates but in general decline faster than the pre-2010 wells. Still, we predict higher EURs for the newly completed wells.
• The more productive newer wells result from recent advancements in completion technology: longer laterals, larger hydrofractures, bigger stimulated reservoir volumes, and more fracture stages.
• Operators have also learned to drill new wells only in the most prolific area of the Bakken region at the center of the Williston Basin.
• With time, negative trends in oil production have amplified in the Bakken. We observe higher GOR values (reservoir degassing); higher water cuts (contacting water-bearing formations and drilling in regions with lower oil saturations); faster decline rates; and excessive well density, especially in the most prolific areas.
• When all most productive areas in the Bakken are extensively drilled, the poor immature areas with high water production will be the only targets left for infill drilling. In this case, technology will not enhance much performance of the future wells.
• We encourage practitioners to adopt our fast and accurate method of predicting oil production in shales that is a viable alternative to the hyperbolic DCA which yields an ‘illusory picture’ of shale oil resources.

Elsewhere [34], we have outlined a hybrid method that combines the precise GEV statistics of horizontal oil wells in the Bakken with the physics-based scaling method developed in this paper and combined with late-time inflow [26]. This hybrid approach follows from our comprehensive evaluation of the Barnett play futures [27].