Gaussian Decline Curve Analysis of Hydraulically Fractured Wells in Shale Plays: Examples from HFTS-1 (Hydraulic Fracture Test Site-1, Midland Basin, West Texas)

Abstract

The present study shows how new Gaussian solutions of the pressure diffusion equation can be applied to model the pressure depletion of reservoirs produced with hydraulically multi-fractured well systems. Three practical application modes are discussed: (1) Gaussian decline curve analysis (DCA), (2) Gaussian pressure-transient analysis (PTA) and (3) Gaussian reservoir models (GRMs). The Gaussian DCA is a new history matching tool for production forecasting, which uses only one matching parameter and therefore is more practical than hyperbolic DCA methods. The Gaussian DCA was compared with the traditional Arps DCA through production analysis of 11 wells in the Wolfcamp Formation at Hydraulic Fracture Test Site-1 (HFTS-1). The hydraulic diffusivity of the reservoir region drained by the well system can be accurately estimated based on Gaussian DCA matches. Next, Gaussian PTA was used to infer the variation in effective fracture half-length of the hydraulic fractures in the HFTS-1 wells. Also included in this study is a brief example of how the full GRM solution can accurately track the fluid flow-paths in a reservoir and predict the consequent production rates of hydraulically fractured well systems. The GRM can model reservoir depletion and the associated well rates for single parent wells as well as for arrays of multiple parent–parent and parent–child wells.

1. Introduction

During the past two decades of intensive field development projects in North American shale acreage (gas and liquids), a plethora of tools and methods has been developed to aid operating companies in the optimization of decisions about well spacing and fracture spacing. One might conclude that the decision-making process is clear by now. However, even with the vast set of tools and methods at our disposal, the production forecasting of shale well systems is still marred by high uncertainty, with a consequent need for improved accuracy of both the well rate predictions and related reserves estimations [1].

While modeling efforts are advancing, honoring the field data by history matching production behavior of wells with the nearly 100-year-old Arps method of decline curve analysis (DCA) [2,3] is still prevalent in our industry as a means of estimating the volume of technically recoverable resources. DCA solutions are fast, cheap and easy to use for computing probabilistic forecasts of future well performance and for the reporting of P90, P50 and P10 resource volumes [4,5].

In this study, a new Gaussian DCA tool, first proposed in Weijermars [6], was applied to analyze the production performance of 11 wells from the Hydraulic Fracture Site-1 (HFTS-1) in the Wolfcamp Formation, a shale play in the Permian Basin, West Texas. The HFTS-1 wells analyzed here were drilled in 2015 and 2016; all were completed in 2016 with zipper-fracking treatment. The treatment details of the HFTS-1 wells were given in prior studies [7,8]. The production performance of the HFTS wells was first analyzed using both the traditional Arps DCA and the new Gaussian DCA tool, and the type curves of the two DCA methods were compared.

The Gaussian DCA method is physics-based and therefore provides certain advantages over Arps hyperbolic DCA (which is currently the most widely used DCA tool in the industry). The Gaussian method can also be used as a forward modeling tool, which is particularly helpful as an alternative or complementary method to both Arps DCA and traditional pressure-transient analysis (PTA).

Three application modes are possible based on new solutions for the pressure diffusion equation using Gaussian pressure transients [9,10]:

• Gaussian DCA is faster than any other DCA method (only one fitting parameter—the hydraulic diffusivity); examples are given in the present study (Section 3).
• Gaussian PTA: After having established the hydraulic diffusivity for a relevant lease domain, one can apply the forward modeling mode and estimate fracture half-length (as in PTA/RTA well-test methods); examples are given in the present study (Section 4).
• Gaussian reservoir models—GRMs work with the full Gaussian solution of the pressure diffusion equation. The Gaussian solution method is closed-form and thus grid-less, arguably more accurate and faster than any of the concurrent modeling platforms used for hydraulically fractured wells. Some examples are given in the present study (Section 5).

Prior studies started out with deriving the Gaussian solution of the diffusivity equation [9] and then developed the solutions to couple the pressure gradient to construct the advancing pressure depletion near hydraulically fractured wells [6]. Next, the theoretical framework was applied to estimate the hydraulic diffusivity for 68 counties with dry gas wells of four major US shale gas plays (Marcellus, Haynesville-Bossier, Utica and Barnett [10]).

Given the unrivalled speed, accuracy and ease of use, the new Gaussian DCA tool and expanded PTA and GRM application modes appear to be very helpful for planning drilling and completion programs for effective field development of shale plays. The capacity for quick iterations and optimizing well design parameters (fracture spacing, well spacing) make Gaussian solutions a mainstay tool. Practical examples are given in the present study (Section 3 to Section 5).

This paper first provides a succinct review of some key aspects of the shale industry as a basis for strategizing whether new tools and methods are still needed. These key aspects are: (1) advances in global shale development (Section 2.1), (2) limitations of key tools and methods used (Section 2.2), (3) key lessons learned over the past decade (Section 2.3) and a review of the HFTS-1 well performance (Section 2.4). Next, the new Gaussian DCA method is briefly explained and applied (Section 3), followed by examples of Gaussian PTA (Section 4) and Gaussian reservoir models (GRM) (Section 5). A discussion follows (Section 6), and conclusions are formulated (Section 7).

2. State of the Art in Shale Basin Development

This paper first briefly reviews four main topics relevant to this paper. These topics are: (1) shale development now taking place in earnest in at least three major hubs outside North America (Section 2.1), (2) limitations of concurrent tools being used for field planning ranging from integrated fracture propagation and reservoir models to traditional tools such as decline curve analysis and pressure-transient analysis (Section 2.2), (3) key lessons learned over the past decade of US shale basin development (Section 2.3) and (4) the performance of wells drilled and studied in the Hydraulic Fracture Site-1 (HFTS-1) in the Wolfcamp Formation, a shale play in the Permian Basin, West Texas, as a prime example of the current state-of-the-art in production analysis and performance forecasting of shale wells (Section 2.4).

2.1. Advances in Global Shale Development

The potentially large hydrocarbon resource volumes harnessed in shale plays around the world was first inventoried in 1997 by Rogner [11]. Shale had already been recognized as a prolific source rock in many US basins, as detailed in local studies such as in the 1997 analysis of gas sorption against kerogen content in the San Juan Basin, Oklahoma [12]. Meanwhile, the development of shale basins in North America (US and Canada) has intensified over the past two decades.

The development of shale basins outside North America initially stalled [13]; among the primary reasons were that, unlike in North America, private landowners under the ruling mineral laws elsewhere are generally not entitled to royalty payments, and densely populated areas are generally not supportive of shale development projects. Meanwhile, decade-long exploration phases of remote shale basins in Argentina, Saudi Arabia and China have been completed; these exploration campaigns are now progressively being converted into massive development projects.

In the Nuequén Basin, Argentina, wells have been completed at a rate of about 200 wells per year since 2014 [14,15], and in the Jafurah Basin, Saudi Arabia, several hundred field delineation wells have been drilled over the past decade, and field depletion will need at least 10,000 wells to be drilled [16,17]. Shale gas production in China started in 2013, with commercial development occurring only in the Wufeng and Long-Maxi formations of the Sichuan Basin, but numerous basins have active exploration programs ongoing [18].

2.2. Limitations of Concurrent Model Solutions

The global drive to open up the development of hydrocarbon resources from shale plays in new regions means that the tools available for such development are continually being evaluated and adapted to the local needs. For unconventional shale plays, the industry has developed fracture propagation tools and production forecasting platforms [19]. The usefulness of advanced tools to forecast the production performance of hydraulically fractured wells remains very limited for three main reasons:

The compounded effect of pragmatic modeling approaches raises concerns regarding the accuracy of the results. These pragmatic simplifications are: (1) violation of the boundary conditions of physical descriptors amalgamated in the computational platform to arrive at seemingly comprehensive outputs, (2) complex parameterization that lacks justification with a sensitivity analysis of those parameters’ importance, (3) upscaling assumptions and grid-refinement limitations leading to approximate results. The compounded effect of these pragmatic modeling approaches on the accuracy of the output remains largely unexplained; no comparative SPE project for testing production forecasting tools of hydraulically fractured well behavior exists.

B.

Non-uniqueness of history matches: Another drawback of commercial modeling platforms is that once a forward model has been constructed, history matching can be achieved by changing a broad number of parameters. It is often claimed that “no one knows what is really happening in the reservoir”. Such simplifications are hard to refute because diagnostic tools are still in their infancy. This also applies to fiber-optics, which provides strain patterns but loses the signal as soon as the next stage is activated; consequently, the state of strain after flow-back remains unexplained [20].

C.

Inability to probabilistically model well behavior: Integrated modeling platforms based on finite grid volumes are computationally demanding and therefore do not lend themselves to probabilistic forecasts where the most critical parameters are inputs based on real-world uncertainty ranges (for example from well-logs). This means even the establishment of P90, P50 and P10 type curves is commonly performed deterministically, using discrete inputs with maximum, minimum and mean values.

The Gaussian solution method of the pressure diffusion equation for hydraulically fractured wells aims to sidestep all the three principal concerns and limitations listed above, as explained and exemplified in Section 3 to Section 5 of the present paper. This claim stems from a modeling approach based on advanced analytical solutions founded in rigorous mathematical descriptions of the time-dependent changes in the reservoir system [6].

2.3. Key Empirical and Theoretical Lessons Learned over the Past Decade

After two decades of intensive drilling in the US unconventional shale basins, one distinct feature has emerged: reduction in the perforation-cluster spacing leads to higher early well rates (Figure 1). The perforation clusters are shot in the cement casing of the wellbore at the start of the fracture treatment with a certain spacing between each cluster; these clusters are the intended points of hydraulic fracture initiation and growth, which may either work or fail (about 25% of the clusters appears to fail; see Waters and Weijermars [21]). However, an important finding from the empirical profiles based on data analytics in Figure 1 is that the effect of fracture-spacing reduction has an immediate impact on the early well rate. Tighter fracture spacing increases the early well rates [22,23].

The relationship between fracture-spacing reduction and well-rate increases has been separately modeled with an analytical flow-cell model, the accuracy of which was benchmarked against independent modeling tools [24,25,26]. All models clearly demonstrated that the closer spacing of the hydraulic fractures does not increase the estimated ultimate recovery but merely leads to accelerated early production (Figure 2a,b). However, the pressure depletion in the reservoir regions between the tighter hydraulic fractures is also accelerated, due to which in such wells the transition to apparent boundary-dominated flow will occur earlier than in wells with wider fracture spacing [25].

Another main finding from prior well model studies was that the reduction in the spacing of wells in shale plays will be reflected in the cumulative production profiles in a unique way; the resulting cumulative production profiles differ distinctly from those resulting from fracture-spacing reduction. Figure 3a,b show such typical cumulative production curves due to well-spacing reduction. Wells with otherwise identical completions but less spacing between them will initially all have the same early well rate. Later in the well life, the wells with the lesser spacing between them will reach true boundary-dominated flow faster than the wider spaced wells. This difference in the advent of pressure interference between wells can be computed: the pressure transient will reach the inter-well boundary earlier in the well life when spacing between the wells is narrower.

Consequently, the effect of fracture-spacing reduction on the well-rate profiles (Figure 1 and Figure 2a,b) can be readily distinguished from the effect of well-spacing reduction (Figure 3a,b). Perforation cluster-spacing reduction (i.e., fracture-spacing reduction) immediately leads to increases in the early well rates (Figure 1 and Figure 2a,b). The other major effect is pressure-depletion interference between adjacent wells due to well-spacing reduction, which is typically seen in production profiles as an effect only later in the well life (Figure 3a,b).

2.4. HFTS-1 Well Performance

The production performance of 11 wells completed in the Wolfcamp Formation at Hydraulic Fracture Site 1 (HFTS-1) is summarized in Figure 4. The HFTS-1 wells were all completed with 660 ft horizontal well-spacing, so differences in well performance cannot be attributed to well-spacing differences. These wells were drilled in 2015 and 2016; all were completed with hydraulic fracture treatment in 2016 using the zipper-fracking mode [7,8]. A new test site (HFTS-2) drilling and research program is ongoing in the nearby Delaware Basin (like the Midland Basin, part of the Permian shale play in West Texas; Zhao et al. [28]), but the HFTS-2 wells were recently drilled and have relatively short production performance; these wells can be analyzed with the Gaussian solution in the future.

What stands out from the HFTS-1 wells is that although they were completed by and large with similar fracture treatment plans, the well performance is highly variable, and the Upper Wolfcamp wells have consistently better production than the Middle Wolfcamp wells (Figure 4). Prior analysis of production and PVT data [7] has revealed that the water cut of the Middle Wolfcamp formation fluid is 3 times higher than that of the Upper Wolfcamp formation fluid, which entirely explains the generally lower oil production rates of the Middle Wolfcamp wells—its production system pumps 3 times more water than oil; in the Upper Wolfcamp, the WOR is 1:1.

The fanning of the cumulative production curves seen in both the Upper and Middle Wolfcamp wells (Figure 4) can be immediately explained as indicative of differences in fracture spacing (and/or equivalent variations in fracture half-length). This can be concluded based on the empirical and theoretical lessons highlighted in Section 2.3. There is no impact of well spacing variations in the case of HFTS-1 wells because all wells have exactly the same inter-well distance (660 ft). The variation in fracture half-lengths that matches the variation in the cumulative production of the HFTS-1 wells was quantified in this study (Section 4).

3. New Method: Gaussian Decline Curve Analysis

The production performance of the HFTS-1 wells was first analyzed here using both the traditional Arps DCA and the new Gaussian DCA tool, and the two methods were compared. Other DCA methods have been proposed during the past two decades of shale development in North America, and these were reviewed and mutually benchmarked in several studies [29,30,31]. Some DCA methods require up to five parameters to be adjusted during the matching process with historic production data, while Arps has three matching parameters (Section 3.1); in comparison, the Gaussian DCA method only requires one matching parameter (Section 3.2).

3.1. Arps Decline Curve Analysis Method

The Arps hyperbolic DCA method is currently the most widely used method for forecasting well production performance based on historic data and uses three variables (q_i, d_i and b) to attain a match:

$$q(t)=q_i\frac{1}{{(1+b{d}_{i}t)}^{\frac{1}{b}}}$$

(1)

The key variables used in Equation (1) are the initial production rate, q_i; the nominal decline rate, d_i and the so-called dimensionless b-value. For the Arps DCA method, one sees that $q(t)$ declines as a function of initial production rate $q_i$ times a hyperbolic function ${(1+b{d}_{i}t)}^{-\left(\frac{1}{b}\right)}$. Although the overall hyperbolic function is dimensionless, it is important to substitute the correct units, with $d_i$ commonly quoted in terms of percentage decline per year but reformulated as a daily or monthly decline rate depending on which the time, $t$, must be correspondingly specified in either days or months. Least-square history-matching on daily data will give decline rates higher than when fitting on monthly data, which is entirely an effect of compounding of the production decline, which appears higher when daily time steps are used.

3.2. Gaussian Decline Curve Analysis

The original pressure in a reservoir is progressively lowered as the pressure transient advances away from a hydraulically fractured well system where a constant bottomhole pressure is assumed. The advance of the pressure transient can be quantified by solving the pressure diffusion equation; a new Gaussian solution method was recently derived from basic principles [6,9]. Gaussian solutions can help to determine the optimum spacing of hydraulic fractures and wells, based on visualizations of the pressure transient advancing from individual hydraulic fractures (Figure 5). The simulator does not require gridding and therefore has unlimited resolution and enables fast computations for optimization iterations.

Gaussian solutions of the pressure diffusion equation also can be simplified into a new decline curve analysis (DCA) method for accurate history matching results on shale well production data. The Gaussian DCA method can match any set of production data, using the following expression [6]:

$$q(t)=q_i\frac{t_1}{t}{e}^{\frac{x^2}{4\alpha }\left(\frac{1}{t_1}-\frac{1}{t}\right)}$$

(2)

Equation (2) specifies the decline of the well rate, $q(t)$, as a function of the initial production rate $q_i$ times the Gaussian decline function $\frac{t_1}{t}{e}^{\frac{x^2}{4\alpha }\left(\frac{1}{t_1}-\frac{1}{t}\right)}$, which features the hydraulic diffusivity $\alpha$ and several scaling parameters such as unit time $t_1$ and unit length scale $x$. Equation (2) was simplified in [6] by introducing two normalizations: $t^{\prime }=t/t_1$ and ${\alpha }^{\prime }=\alpha t_1/x^2$. Substitution of these normalized parameters into Equation (2) gives:

$$q(t^{\prime })=q_i\frac{1}{t^{\prime }}{e}^{\frac{1}{4{\alpha }^{\prime }}\left(1-\frac{1}{t^{\prime }}\right)}$$

(3)

Dropping the asterisks from Equation (3), we obtain Equation (30) of Weijermars [6]:

$$q(t)=q_i\frac{1}{t}{e}^{\frac{1}{4\alpha }\left(1-\frac{1}{t}\right)}$$

(4)

The Gaussian DCA function of Equation (4) can history-match real well production rate data using only one variable (the hydraulic diffusivity, $\alpha$) for the history matching. The time, t, initiates automatically when first production starts; it is therefore not an arbitrary or independent system variable but simply starts from zero. Likewise, the initial production rate, q_i, will in the physics-based Gaussian DCA model also start from nearly zero (a practical initial unit volume rate is used in our models as q_i is the rate per time unit $t_1$); the actual initial well rate starts would start from zero, because no pressure gradient exists until the pressure transient has managed to travel into the reservoir. However, at the moment the pressure transient starts to move into the reservoir, which is when the pressure gradient comes into physical existence and Equation (4) applies, which is why using a small, initial rate of unit production is warranted.

In summary, the only unknown key variable in Equation (4) is the hydraulic diffusivity, $\alpha$, which is scaled in numerical values determined by the chosen units for time, $t_1$, and units for length, $x$ [used in the normalization of Equation (2)]. For example, if the assumed time unit $t_1$ is in unit days, and we history match oil well production rates, then the production rate numbers used for q is in bbls/day with q_i = 1 bbl/day at time zero. Likewise, if we simply settle on a length scale $x$ of unit ft, then the units of α will be in ft²/day (in accordance with the time unit adopted for $t_1$). In the present study, $t_1$ = 1 day and $x$ = 1 ft; which is recommended for all Gaussian DCA analysis. For gas wells, q is in Mcf/day with q_i = 1 Mcf/day at time zero [10].

Having just one history matching parameter ($\alpha$) in the Gaussian DCA method of Equation (4) gives excellent history-matching results to the historic well data, and isfaster than the Arps DCA method of Equation (1), which requires convergence of the matching DCA curve and the real data by a least-square fit using three parameters (q_i, d_i and b). Gaussian DCA also uses least square fitting, but only one parameter ($\alpha$) needs to be solved. Examples of such matching results are given in the next section.

3.3. Matches for HFTS-1 Well Rates with Gaussian DCA and Arps DCA

Well-rate versus time history-matches were obtained for all 11 HFTS wells. Typical matches with both the Gaussian DCA and Arps DCA methods are shown for Well 6U in Figure 6a,b. The Gaussian curve matches the early production data very closely (Figure 6a). However, none of the prior DCA methods can match early well data. This is also the case with the Arps DCA method (Figure 6b) because it is not based on a valid solution of the pressure transient (unlike the Gaussian DCA method); therefore, Arps-based decline curves never start producing the well from a zero rate but instead start out as an asymptote coming from infinity or a large initial well rate. This is physically incorrect—all wells only come on stream slowly. The Gaussian DCA closely matches the waxing of early well rates until a peak rate is reached (Figure 6a), after which the well rate declines asymptotically.

The Gaussian DCA solution is based on well physics because it was obtained by solving the pressure diffusion equation for the pressure transient using basic principles [6]. Consequently, the Gaussian DCA method is not only faster but also more accurate than the Arps DCA method in that the Gaussian method accurately shows the onset of early production from a zero rate at time zero (Figure 6a); at that moment the pressure-transient does not exist because prior to first production the bottom hole pressure in the well system is not lowered yet. In contrast, the Arps hyperbolic DCA method starts at an arbitrary large production rate at time zero (Figure 6b), because it is a trend line, with little or no connection to the well physics.

Appendix A gives history-matched curves for all of the HFTS-1 wells, giving side-by-side results using both the Arps and Gaussian DCA methods. The corresponding DCA matching parameters are also summarized in Appendix A. Applying the Gaussian method to real production data of the HFTS-1 wells affirmed the practicality and accuracy of the new Gaussian DCA method.

4. Gaussian Pressure-Transient Analysis (PTA)

After having obtained the hydraulic diffusivity for a relevant lease domain (as summarized for HFTS-1 wells in Appendix A,

Table A2. Gaussian DCA matching parameters for Upper (left) and Middle (rights) Wolfcamp wells.

Well	bbls/Day	ft2/Day	m2/s	Well	bbls/Day	ft2/Day	m2/s
	qi	α	α		qi	α	α
3U	1	0.0222	2.38 × 10−8	4M	1	0.0225	2.42 × 10−8
4U	1	0.021777	2.34 × 10−8	5M	1	0.022487	2.42 × 10−8
5U	1	0.022035	2.37 × 10−8	6M	1	0.022843	2.46 × 10−8
6U	1	0.022243	2.39 × 10−8	7M	1	0.022243	2.39 × 10−8
7U	1	0.0223	2.39 × 10−8	8M	1	0.0224	2.41 × 10−8
8U	1	0.0225	2.42 × 10−8

), one can apply a fuller modeling mode based on Gaussian solutions of the pressure diffusion equation [6]. The simplest mode of use was Gaussian DCA (see the three modes cited in the Introduction of this paper), but now we advance to history matching well data with the Gaussian pressure-transient solution for Gaussian PTA, which can be used to estimate fracture half-lengths (as in the traditional PTA/RTA approaches used in well-test methods). Examples are given below. First, the theory is briefly explained (Section 4.1), then results are given (Section 4.2).

4.1. Key Equations of Gaussian PTA

The total influx of reservoir fluid, $q_r(t)$, with viscosity, $\mu$, from a reservoir space with permeability, $k$, and porosity, $\varphi$, into the collective system of hydraulic fractures with total surface area, $A$, (flux from two sides gives the term $2A$) is ([6], Equation (28)):

$$q_r(t)=2A\varphi \frac{k}{\mu }\frac{P_0-P_{BH}}{2t\alpha }x{e}^{-\left(\frac{x^2}{4t\alpha }\right)}$$

(5)

Equation (5) further features as inputs: the hydraulic diffusivity, $\alpha$, and then shows how the well rate will decline over time, $t$, for a well system produced with a pressure differential of $P_0-P_{BH}$ (with bottomhole pressure $P_{BH}$ for the well system and initial pressure $P_0$ for the unperturbed reservoir). The length scale is again conveniently fixed by the scaling parameter, $x$, which must be assigned the same value as used in the DCA normalization for estimating $\alpha$(when using field units, $x$ = 1 ft is recommended).

Equation (5) gives the total fluid rate in reservoir barrels, which translates to the rate of stock tank barrels of total fluid, $q_w(t)$, produced at the wellhead as follows:

$$q_w(t)=\frac{q_r(t)}{B_0}$$

(6)

with the formation volume factor $B_0$ (bbl/stb) accounting for the loss of volume at the surface due to the depletion of gas saturation; correspondingly, $B_g$ (Mcf/Mscf) can be used for Gaussian RTA on gas wells.

4.2. Gaussian PTA Results

Having solved the hydraulic diffusivity, $\alpha$, with the Gaussian DCA method for Well 6U (Figure 6a, by history matching using Equation (4)), one may now infer the fracture-half-length that can sustain the observed well rate by again history matching the production rate of Well 6U by applying Gaussian PTA using Equations (5) and (6). The total fracture area, $A$, is estimated first and is related to average fracture half-length, W_f, (assuming all n perforation clusters are successful and produce n hydraulic fractures) and the pay-zone height, h, as follows:

$$W_f=\frac{A}{2nh}$$

(7)

The area, $A$, used in our models is assumed to correspond to the surface area of only those fractures which significantly contribute to the well rate by having effective hydraulic conductivity close to infinity.

Table 1. Input parameters used in the Gaussian PTA history-matching.

Reservoir Attribute	Unit	Value
Initial Pressure UWF	psi	4073
Initial Pressure MWF	psi	4250
Bottomhole pressure	psi	500
Porosity	none	0.06
Permeability	Darcy	5.0 × 10−7
Viscosity	cPoise	0.6
Payzone height	ft	100
Formation volume factor	bbl/stb	1.4186

lists the required input parameters used in the Gaussian RTA history-matching of the daily production performance of all HFTS-1 wells. There is no uncertainty in the values of $P_0$ and $P_{BH}$. For the wells studied here, the porosity of the reservoir rock can be estimated from well logs, which likewise involves little uncertainty. The viscosity of the fluid system under reservoir conditions and the formation volume factor above the bubble point can be inferred from PVT reports, which therefore also result in negligible uncertainty regarding what these input values should be.

The principal uncertainty remaining for the output of Equation (5) is the permeability, $k$, of the stimulated rock volume. A deterministic approach, using a discrete permeability value of 500 nanoDarcy, results in estimations of the fracture half-length, W_f, for the Upper and Middle Wolfcamp wells, as summarized in

Table 2. Result for fracture half-lengths of Upper Wolfcamp wells from Gaussian PTA history-matching.

Unit	Well Number→	8U	7U	6U	5U	4U	3U
Darcy	Permeability	5.0 × 10−7	5. 0 × 10-−7	5.0 × 10−7	5.0 × 10−7	5.0 × 10−7	5.0 × 10−7
ft2/day	Diffusivity	0.0225	0.0223	0.022243	0.022035	0.021777	0.0222
m2/s	Diffusivity	2.42 × 10−8	2.39 × 10−8	2.39 × 10−8	2.37 × 10−8	2.34 × 10−8	2.38 × 10−8
ft	Height	100	100	100	100	100	100
ft	2Wf	816	788	1062	1167	802	1102
ft	Wf	408	394	531	583	401	551
Clusters	Fractures	113	149	113	113	186	113

and

Table 3. Result for fracture half-lengths of Middle Wolfcamp wells from Gaussian PTA history-matching.

Unit	Well Number→	8M	7M	6M	5M	4M
Darcy	Permeability	5. 0 × 10−7	5.0 × 10−7	5.0 × 10−7	5.0 × 10−7	5.0 × 10−7
ft2/day	Diffusivity	0.0224	0.022243	0.022843	0.022487	0.0225
m2/s	Diffusivity	2.41 × 10−8	2.39 × 10−8	2.46 × 10−8	2.42 × 10−8	2.42 × 10−8
ft	Height	100	100	100	100	100
ft	2Wf	947	625	776	900	557
ft	Wf	473	313	388	450	279
Clusters	Fractures	113	183	113	113	185

, respectively.

5. Gaussian Reservoir Models (GRMs)

Normalization of a new Gaussian solution of the pressure diffusion equation gave the Gaussian DCA formula [Equation (4)], which is useful for determining the hydraulic diffusivity responsible for those well rates in history matches with actual well data (Section 3). Once the hydraulic diffusivity is established for the rock volume drained by the well system, it is possible to apply Gaussian PTA to establish the fracture half-length of each HFTS-1 well (Section 4). Next, Gaussian pressure transients can be applied to evaluate the production rate (including pressure-interference effects)for for hydraulically fractured well systems prior to drilling with the full set of Gaussian solutions (GRM mode of modeling; see the three Gaussian application modes in Section 1). Figure 7 shows an example of a comprehensive model with (1) pressure depletion, (2) fluid velocities and (3) flow paths in the reservoir region that is drained by the well system. Two key solutions, one approximate and one exact, derived in the prior companion study [6] are summarized and commented on below to support their practical use in field applications. These Gaussian solutions and applications are not limited to shale wells but are equally valid for geothermal well systems (in which case thermal diffusivity solutions given in [6] can be included).

Probabilistic Reserves Estimations

For shale operators it is of paramount importance to comply with reserves estimation guidelines, which requires the determination of proved reserves. Increasingly, probabilistic methods are used to unequivocally establish P90 proved reserves to distinguish from unproved reserves (comprising the probable (P50–P90) and possible (P10–P50) reserves categories), which in the US may not be used for reserves-based capital lending to finance the cost of new wells [32]. Distinguishing between P90, P50 and P10 reserves estimations using deterministic inputs often leads to ambiguity whether ‘best’ estimate indeed is a P90 value, which is why the uncertainty range in the reserves can be more appropriately constrained by applying probabilistic methods.

Providing probabilistic reserves estimations (for a given day of first production forward) for new wells to be drilled in the same acreage is possible with Gaussian methods, applying two possible modeling modes:

(1)

Using Gaussian DCA in forward modeling mode with probabilistic inputs. For example, a probability density input function can be created from the hydraulic diffusivity-values ($\alpha$) of Table A2.

(2)

Using Gaussian PTA in forward modeling mode with probabilistic inputs for the fracture half-lengths as well as the probabilistic inputs for $\alpha$. A probability density input function can be created from the fracture half-lengths estimated for the existing wells in the acreage (for example, as per Table 2 and Table 3 in Section 4.2).

The input distribution for $\alpha$, used for reserves estimation according to the first method, is given in Figure 8a. The corresponding output of the probabilistic reserves volumes is charted in Figure 8b. The cumulative production curves for P10, P50 and P90 wells are given in Figure 8c. Of course, the well performance forecasts and estimated ultimate recovery (EUR) estimations assume that the new wells will be completed with the same technology and fracture spacing as was used in the existing wells based on which the estimations were made. If an operator chooses to apply different fracture spacing and different fracture half-lengths (due to technology innovation), then these inputs suffice to construct a forward model with Gaussian PTA (using the Equations of Section 4) that can still be made probabilistic by using only a probability distribution for the history-matched hydraulic diffusivities; the fracture dimensions are then fixed according to the engineering design for the hydraulic fractures.

6. Discussion

The need for optimization of the well performance of hydraulically fractured wells in shale plays has led to the development of a vast range of tools and methods to support professionals in operating companies. It was argued that concurrent model solutions still suffer from considerable inaccuracy due to a variety of causes. For example, DCA methods are not physics-based and/or violate initially assumed boundary conditions of boundary domination flow (as is the case in Arps method applied to shale wells). Concurrent commercial simulators are all grid-based, which (among other drawbacks) makes them ill-suited for a probabilistic forecast of well-behavior (Section 2.2). The new Gaussian solution methods presented here open new avenues for practical applications. Some important directions and insights are highlighted in the five discussion items of the following subsections.

6.1. Workflow Recommendations

The set of tools presented here can be readily used in a workflow where each previous step contributes to reduce uncertainty of input parameters for the next step. The flow chart of Figure 9 summarizes the proposed workflow. Early steps in the workflow are aimed at reducing the uncertainty in unknown parameters by history matching; the history-matched parameters can then next be used in forward models.

The Gaussian DCA tool mode is the least demanding on data inputs; just the well rates (using daily rates rather than monthly rates is recommended) are sufficient for completing Step 1. Once the hydraulic diffusivities for a bunch of wells are estimated, they can be used as probability distribution input for hydraulic diffusivity (see example in Section 6.1) to produce probabilistic reserves estimations in Step 2. Switching to the Gaussian PTA tool mode, Step 3 can be used to constrain—for each well—the average fracture half-length. Once those half-lengths are found, they can be used as probability distribution inputs in addition to the probabilistic inputs for hydraulic diffusivity to produce probabilistic reserves estimations in Step 4.

The reserves estimations of Steps 2 and 4 need not be the same, but as long as the prescribed workflow can be executed by different engineers leading to the same outcomes in each step, then the reserves reporting using the Gaussian methods (Step 2, Step 4 or combinations) will be valid under existing reserves reporting guidelines. Step 5 is a more elaborate approach where the effect of the production and well system placement on the reservoir can be systematically investigated. Of course, coupling with discounted cash flow models is possible in all steps to compute which fracture spacing and well placement yields the highest return on investment.

6.2. Strength and Weaknesses

Gaussian solution methods are a new and wholly original avenue for solving well performance and reservoir management issues. The term Gaussian may sound familiar, but no prior solution of the pressure diffusion equation existed without reference to any well rate as an input parameter in the well-testing literature. The solution of the pressure-transient without reference to the well rate as an input was first given in Weijermars [6,9]; this independence of well rate is one of the major strengths and innovations of the new solutions for the pressure-transient. The present article provides practical examples of how the new tools can be applied in practical applications. As of the date of this study, no obvious shortcomings of the methods have become apparent, which is why the suggested methodologies are assumed to become mainstay staple tools for petroleum engineering operations.

6.3. Future Work

The new methods (DCA, PTA. GRM) for reservoir and well performance analysis presented in this study, based on new solutions of the pressure-transient equation, can be applied in a vast range of subsurface reservoir utilization projects: water disposal, hydrogen storage, carbon-dioxide sequestration, geothermal extraction and hydrocarbon extraction projects (not being limited to shale plays only). To further pave the way for such utilization, more case study examples are planned for future work by the author and his research team. Separately, a company may be formed with the aim of making available a practical suite of tools for easy usage in practical operations based on the various applications of the Gaussian solutions (see Figure 9), of which some examples were provided in this paper.

7. Conclusions

This study offered practical examples of how the key Gaussian equations based on rigorous solutions of the pressure diffusion equation in the porous medium of a subsurface reservoir, produced with hydraulically fractured wells, can be used in practical applications. The new solutions for the Gaussian pressure transients also yielded new formulae for physics-based Gaussian decline curve analysis (DCA), suitable for one-parameter production history matching, production forecasting and reserves (EUR) estimations. The reserves estimations can be based on probabilistic methods, which gives accurate estimations for P10, P50 and P90 EUR volumes. All estimations are very accurate because they are based on closed-form solutions of the pressure diffusion equation. This study also includes practical examples of Gaussian pressure-transient analysis (PTA) for estimating fracture half-lengths. Hydraulic diffusivity estimations are given for the drained reservoir volume. Full vector field solutions for pressure, velocity and flow paths are possible with the Gaussian reservoir modeling mode.