The Role of Natural Fractures in Shale Gas Production

Natural fractures seem to be ubiquitous in shale gas plays. It is often said that their presence is one of the most critical factors in defining an economic or prospective shale gas play. Many investigators have presumed that open natural fractures are critical to gas production from deeper plays such as the Barnett, as they are for shallower gas shales such as the Devonian shales of the northeastern US and for coal bed methane plays. A common view on production mechanisms in shales is “because the formations are so tight gas can be produced only when extensive networks of natural fractures exist” [6]. However, there is now a growing body of evidence that any natural fractures that do exist may well be filled with calcite or other minerals and it has even been suggested that open natural fractures would in fact be detrimental to Barnett shale gas production [9].

Commercial exploitation of low mobility gas reservoirs has been improved with multi-stage hydraulic fracturing of long horizontal wells. Favorable results have been associated with large fracture surface area in contact with the shale matrix and it is here that the role of natural fractures is assumed to be critical. For largely economic reasons hydraulic fracturing for increasing production from shale gas reservoirs is often carried out using large volumes of slickwater injected at pressures/rates high enough to create and propagate extensive hydraulic fracture systems. The fracture systems are often complex, due essentially to intersection of the hydraulic fractures with the natural fracture network. After hydraulic fracturing operations the injected water is flowed back. Typically, only a small percentage (on the order of 20 to 40%) is recovered.

In this paper we investigate the role played by natural fractures in the gas production process. By applying a new model of the production process to data from many shale gas wells across a number of shale plays in North America, we can for the first time begin to sort out assertion from inference in the role that these fractures play. Specifically, we are able to estimate the magnitude of the fracture surface through which gas is actually produced. We are able to demonstrate that although it may be commensurate with the expected surface area of open natural fractures for the ultra-low permeability shallow gas shales, it is in fact commensurate with a very much smaller area for the deeper gas shales such as the Barnett. Furthermore, given a typical value of the matrix permeability, almost all the gas between the fractures would have been produced in an uncharacteristically short period of time unless the producing fractures are 100s of feet apart. The implications of these findings for completion and stimulation strategies will be discussed.

1. Introduction

One of the main premises of our investigation of the production processes in shale gas plays is that the industry’s mental picture of the process remains very much influenced by the concepts developed in the 1990s of the production processes in coal bed methane resources and in shallow shale gas plays. It is appropriate therefore that we begin our discussion of shale gas production characteristics by reviewing this early work.

1.1. Coal bed methane

Coal is a heterogeneous and anisotropic porous medium, characterized by two distinct porosity systems:

• Micropores of diameter of the order of nanometers with almost zero permeability.

• Macropores or cleats, slot-like with spacing of the order of 2 cm and width of the order of microns; permeability is stress-dependent but far in excess of the micro-pore permeability. They are often formed by shrinkage of the coal matrix due to dewatering during the coalification process.

Gas is stored essentially by adsorption in the coal matrix; very little is stored as free gas in the micropores or as free gas or dissolved gas in the connate brine in the cleats. In the subsurface the cleats are usually filled with water, some of which must be produced to the surface to facilitate gas production.

The conventional view of the production process divides it into three stages (see Figure 1):

1. Cleat dewatering, lasting of the order of several years; the water production rate gradually falls as water is removed from the cleats. At the same time more and more gas is produced at increasing rates and the relative permeability to gas in the cleats increases leading to lower pressures and more gas production.

2. Stabilized flow: eventually, most of the water in the cleats has been removed, the cleat fluid pressure bottoms out and the relative permeability to gas levels off. Over this period the gas rate slowly peaks.

3. Decline: there is then no more increase in drawdown available to sustain gas production and gas production declines. If the cleat pressure was constant and the pseudo-steady-state (PSS) regime is applicable (see the later discussion for definitions of this flow regime), then gas production rate should decline exponentially.

Figure 1.

Stages in gas and water production from coal (after [1]).

There are three essential elements to a model of the CBM production process:

1. Transport in the coal matrix, modeled as a diffusive process using Fick’s diffusion law. In principle the gas concentration in the coal matrix satisfies a diffusion equation, but it is common to use a pseudo-steady-state (PSS) approximation similar to that proposed by Warren and Root [2] in their dual porosity formulation of production from naturally fractured reservoirs. For example, King et al [3] used the PSS simplification to reduce computing time and because after a period of time the numerical accuracy was deemed to be quite acceptable. We have estimated from King’s data that the PSS solution is valid beyond about 40 days, which is much shorter than the typical duration of the production process. We note, however, that this time scale depends on the assumed values of the diffusivity in the matrix and on the spacing of the cleats, assumed to be of the order of a few cm.

2. Desorption at the cleat/matrix interface as characterized by the Langmuir isotherm

3. Transport of water and free gas in the cleat system. To avoid difficulties in defining the configuration of the cleat system, it is common to adopt a dual porosity description in which the cleat system is treated as a continuum with system characteristics analogous to those of a porous medium. Two-phase flow in this system can for the most part be adequately described by Darcy flow. In narrower cleats it may be necessary to include capillary pressure and slippage effects especially at low pressure.

It is apparent that the cleat or natural fracture system plays a very important part in the production process. The density of the cleats plays two critical roles: first, the close spacing of the cleats reduces the time required for the gas to diffuse to the cleats and, second, it is associated with a high cleat/matrix surface area without which economical gas production would be unlikely. The width of the cleat is a primary influence on the pressure drop in the whole system and therefore on the water and gas production rates. In situ the cleats are usually water-filled and presumably kept open by the pressure of the fluid they contain. The cleats may close somewhat as the pressure falls during production, though this may be more than offset by matrix shrinkage as the gas is desorbed.

1.2. Devonian shales of the appalachian basin

Gas production from Devonian shales received a great deal of attention in the 1980s and early 1990s as a result of US DoE initiatives. This is well documented in many GRI reports and industry publications. The consensus view is that these reservoirs are highly fractured containing a substantial number of fractures with spacing of the order of 1-10 cm (see, for example, [4]). Luffel et al [5] measured the matrix permeability at less than 0.1 nd. Water content of the Devonian Shale averages 2.5 to 3% of bulk volume and appears to be at irreducible water saturation. Typical depth is a few thousand feet, pore pressure is less than about 3000 psi and about 50% of the gas in place is adsorbed; there is little or no water production.

Carlson and Mercer [6] summarized the consensus view of the production process as “because the formations are so tight gas can be produced only when extensive networks of natural fractures exist.” The extent to which this statement holds for other gas shale plays is debatable, but it has certainly been influential in developing the industry’s vision of what is happening downhole.

Gatens et al [7] used a dual-porosity model similar to that formulated by Warren and Root [2] but extended to use the unsteady-state equation instead of the pseudo-steady-state (PSS) equation for matrix flow. Analysis of hundreds of Devonian wells showed that most of the production data fell into the linear transient regime (as we discuss later in this document). Luffel et al [5] obtained a good history match with data by assuming an open fracture spacing of a few feet, while Carlson and Mercer [6] needed a fracture spacing of about 80 ft, both with matrix permeability of less than 0.1 nd. An issue that does not seem to have been addressed however, is whether these fractures, if present, are in fact open and if so how they are maintained open against closure stress. If, like coal cleats, they are initially water-filled, is water production observed? An implicit assumption seems to be that they are open and gas-filled.

Carlson and Mercer [6] proposed that molecular diffusion is the dominant transport mechanism in the matrix in these extremely tight reservoirs, in which case a matrix diffusion coefficient should be used instead of the matrix permeability. They did not evaluate the consequences of this hypothesis. It remains a possibility that the use of such a coefficient would reduce the need for a large fracture surface area and ultimately for the need to propose the existence of a large open fracture network.

Thus there are three essential elements to a model of the production process in the Devonian shales:

1. Desorption of gas in the matrix (as characterized typically by the Langmuir isotherm)

2. Transport in the matrix towards the fracture network, modeled as Darcy flow even though the permeability is extremely small.

3. Transport of free gas in the fracture system.

1.3. Devonian Antrim Shale of the Michigan basin

The Antrim Shale is a shallow, under-pressured, naturally-fractured shale reservoir with characteristically low matrix permeability, and with adsorbed gas, free gas and mobile water co-existing in the reservoir. A typical Antrim well will produce considerable quantities of water early in life, and as dewatering of the reservoir progresses, water production rates decline and a corresponding increase in gas production is normally observed (as a result of gas desorption with reduced reservoir pressures), similar to a CBM well. In fact, the Antrim shale is often considered to be a hybrid of productive dry gas shale and CBM plays. It has characteristics which are similar to these other unconventional reservoirs, but it is also different in many ways. The Antrim shale is more intensely fractured than the Devonian Shales of the Appalachian Basin, with fracture spacing as close as 1 to 2 ft. Kuuskraa et al [8] have noted that the “intensity and interconnection of the fractures govern the shale’s natural producibility.”

The typical depth of the Antrim shale is less than about 2000 feet, pore pressure is a few hundred psi and more than 70% of the gas in place is adsorbed, the remainder being stored essentially as free gas in the matrix pores. Peak gas may occur as late as 3 years into production. Production data has been history-matched using similar software to that used for CBM [8]. It was found that fracture spacing of the order of a few inches facilitated a good match with production data. It was stated that if a fracture spacing of 3-6 ft was used (which is compatible with observations from cores and logs), then production would be an order of magnitude lower than observed in existing wells. One possible resolution of this conundrum may lie, as the authors suggest, in detrital silt layers within the matrix that could provide conductive flow paths. (An alternative explanation that remained unconsidered by the authors lies in the use of the PSS approximation for matrix transport, which may be completely invalid in this context.)

This description leaves many issues unresolved, but importantly places the Antrim shale as an intermediate between CBM and the other shales, in that the natural fractures appear to be conductive and initially water-filled, but has free pore gas in addition to adsorbed gas.

1.4. Deeper gas shales

Most modern commercial gas shale plays are similar in many respects to the Barnett shale, though there are of course many differences and variations in the values of the parameters that control the gas production process. In these relatively deep and high pressure reservoirs, most of the gas is stored as pore gas, but the production process is still similar to that described above for shallower shale gas plays.

The subject of open natural fractures is one of the most contentious within the community of Barnett workers. Many investigators have presumed that open natural fractures are critical to Barnett gas production, as they are for the shallower gas shales, even though there is now a growing body of evidence that any natural fractures that do exist may well be filled with calcite or other minerals (see Figure 2). There are also arguments that suggest that if there was an abundance of open natural fractures within the Barnett, there would be a much smaller gas accumulation present within the reservoir. Open natural fractures, if they existed, would have led to major expulsion and migration of gas out of the shale into overlying rocks, substantially decreasing pore pressure within the Barnett and, hence, the amount of gas in place. The Barnett would not be over-pressured (that is, over-pressured relative to the bounding strata) if copious open natural fractures existed. Note that the Barnett is not just the gas reservoir, but also the source, trap, and seal for the gas; if the seal is fractured and inefficient, then the present gas in place would be reduced because the free gas would be lost, and only the adsorbed gas would remain in the shale (a similar situation to that of the Antrim Shale of northern Michigan). The huge amount of gas in place, in an over-pressured and fully-saturated (in terms of sorption) state, is ultimately what makes the Barnett so prolific.

Figure 2.

Mineralized natural fractures in a Barnett shale sample (adapted from [9]).

A common argument for the necessity for open natural fractures in shale gas plays is that a large surface area is necessary for economic gas production from these very tight rocks. Later in this paper we analyze production data to estimate the magnitude of the fracture surface through which gas is actually produced. We are able to demonstrate that although it may be commensurate with the expected surface area of open natural fractures for the ultra-low permeability shallow gas shales, it is in fact commensurate with a very much smaller area for the deeper gas shales such as the Barnett. Furthermore, given the typical permeability of the Barnett shale (some 100 times that of the shallow gas shales), almost all the gas between the fractures would have been produced in an uncharacteristically short period of time unless the fractures are 100s of feet apart. These issues and conclusions will be discussed at length later in the paper.

2. Production mechanisms, production modeling techniques and simulators

Having outlined the pertinent characteristics of unconventional gas reservoirs, we now document the likely production mechanisms in the various shale plays based on our understanding of their geology and the underlying geophysics.

The matrix permeability of shale gas reservoirs is extremely small, probably on the order of one tenth of a microdarcy or 100nd. It is virtually impossible to produce gas from these reservoirs in commercial quantities unless the wells are hydraulically fractured and even then, or so it is commonly believed, production is really only possible because a network of natural fractures is opened up. (It is interesting to note that gas has been produced from the ultra-tight Devonian shale plays of the North Eastern USA from more conventionally-fractured vertical wells, which implies that multi-stage hydraulic fracturing was unnecessary for these plays. This is the first hint that the role of the natural fractures may be quite different for the Devonian plays and the deeper shale plays.)

An essential element of a mathematical model of gas production from shales is therefore the ability to describe flow in a very tight rock matrix and flow in a network of fractures. In most gas shale reservoirs most of the reservoir fluid is stored in the matrix and the primary flow path is from the matrix into the fractures and thence into the wellbore. There are essentially two methods of characterizing a multiply-fractured reservoir:

• Discrete fracture network (DFN) model, in which the fractures are defined explicitly in terms of their location in the reservoir, their connectivity to one another and to the wellbore and their production characteristics, such as permeability and conductivity.

• Dual porosity/dual permeability models in which the fracture network is treated as a continuum in much the same way as a porous medium is treated as a continuum for analysis of flow characteristics.

2.1. Discrete fracture network models

Many commercial numerical reservoir simulators have the capability of simulating flow through a complex network consisting of pores and fractures. However, one of the greatest drawbacks and limitations of simulating a discrete fracture network model is the a priori assumption that all relevant properties of the fracture network are known. Nevertheless great insights can be obtained into the impact of the essential physical processes by examining simple fracture configurations. We note that in principle many different physical and petro-physical components can be included in numerical simulations. However, in practice it is quite common to see results presented only for the special cases:

• Reservoir fluid of small and constant compressibility.

• Production under constant drawdown conditions.

• No desorption.

• Darcy flow in fractures and matrix.

• Matrix and fracture permeability independent of pressure; it is often assumed that fracture conductivity is essentially infinite.

The simplest fracture network that has been applied to shale gas production consists of a number of planar fractures placed transversely to a horizontal wellbore as illustrated schematically in Figure 3. It is apparent from many published numerical studies that under these circumstances flow from the reservoir can be described in terms of a number of identifiable flow regimes. The following account is taken from a recent paper by Luo et al [10]. These authors used a commercial reservoir simulator to calculate the flow into a horizontal well with six infinitely-conductive transverse fractures as shown in Figure 4.

Figure 3.

Idealized discrete fracture network showing multiple transverse fractures originating from a horizontal well.

Figure 4.

Numerical solutions for flow into six infinitely-conductive transverse fractures (taken [10]).

The flow behavior can be conveniently discussed in terms of five flow regimes as follows.

• Bilinear or linear flow: soon after the well is placed on production reservoir fluid flows normal to the fracture planes and along the fractures into the well. The streamlines are shown in yellow in Figure 4; reservoir pressure is in red and the constant bottomhole or well pressure is in blue. Note that flow into the fracture tips is negligible and each fracture behaves independently of the other fractures. This regime may also be termed the infinite-acting regime in the sense that the neighboring fractures are effectively at infinity. The duration of this regime depends, as we shall see later, on many parameters including the matrix permeability, the fluid compressibility and the fracture spacing. The illustrations in Figure 4 are for infinite conductivity fractures.

• Early radial/elliptical flow: flow into the fracture tips is present, but weak; flow into the fracture surfaces is still predominantly linear, but fracture interference is just beginning to impact the flow. Note that the pressure drawdown in the matrix has almost reached the mid-line between the fractures. At this point the flow regime may be described as pseudo-steady-state or fracture-boundary-dominated.

• Compound formation linear flow (CFL flow): here the fractures are fully interactive and the reservoir drainage area is dominated by the area defined by the length of the well and the length of the fractures. Flow from beyond this area grows in importance.

• Pseudo-radial/elliptical flow: flow from beyond the wellbore/fracture area grows in importance and appears to be radial or elliptical.

• Reservoir-boundary-dominated flow: ultimately the outer boundary of the reservoir begins to impact the flow.

It is difficult to infer from these simulations the time scales and duration of these flow regimes for parameter values other than those used in the particular simulation. Indeed, this highlights one of the severe drawbacks to the full numerical approach to modeling production flow in these reservoirs: it is difficult to make general conclusions about the characteristics of the flow and their dependence on the input data without undertaking very many numerical simulations; this is a formidable task even for a restricted input data space. However, based on the semi-analytic models that are described below, we believe that for many shale gas wells it would be unusual to expect to encounter the compound formation linear flow regime for at least 10 years after the well was placed on production.

2.2. Dual porosity/dual permeability models

The conventional view of a naturally-fractured reservoir is that it is a complex system composed of irregular matrix blocks surrounded by a network of more highly permeable fractures. In reality in tight gas shales some or most of the fractures may not be open to flow or they are opened up only during the hydraulic fracturing process. Warren and Root [2] were among the first to recognize that the simple model of reservoir flow based on single values of the permeability and porosity does not apply to naturally-fractured reservoirs, though they had in mind reservoirs quite different from gas shale reservoirs. In order to handle the problems associated with lack of detailed information on the structure of the fracture network they proposed a dual-porosity model in which a primary porosity associated with inter-granular pore spaces is augmented by a secondary porosity related to that of the network of natural fractures. At each point in space there are two overlapping continua—one for the matrix and one for the fracture network. The detailed geometry of the fracture system need not be specified in this model, but can include as particular examples any of the discrete fracture models described above. In typical shale gas applications the matrix has high storage capacity but low permeability and the fractures have relatively low storage capacity and higher permeability. It is quite possible (or, indeed, likely) that in many shale gas reservoirs no gas is stored in the fractures, though they may become filled with frac fluid during the hydraulic fracturing process.

In the dual porosity formulation flow from the matrix to the fractures is described by a transfer function with Darcy flow characteristics. The original Warren-Root models incorporated the pseudo-steady-state assumption in the matrix blocks and assigned a single value to the pressure within the blocks; the mass transfer rate from the matrix to the fractures depends then on the pressure differential between the matrix and the fracture. Thus these models assumed, almost implicitly, that sufficient time had elapsed that the flow in the matrix blocks between the fractures was already fracture-boundary-dominated. Later in this paper we estimate the time scale on which inter-fracture pseudo-steady-state begins and find that it is typically of the order of several years for a fracture spacing of 100 ft or more. This is quite consistent with typical simulation results described above. If the fracture spacing was as small as 10ft, we should expect to see fracture interference or the onset of PSS flow after about 10 days. For shale gas reservoirs, more complex models (unsteady state or fully transient) are needed to resolve the flow in the matrix in more detail.

A detailed discussion of the Warren-Root model, its background and similar contemporary models can be found in the excellent monograph by van Golf-Racht [11]. We note in particular that Kazemi [12] was one of the first to extend the Warren-Root model to include transient flow in the matrix blocks. Some seventeen years after Warren and Root published their seminal work, Kucuk and Sawyer [13] adapted their model for flow in shale reservoirs by incorporating effects such as desorption from the organic matrix material and Knudsen flow in the pores and, of course, incorporating full transient effects in the matrix blocks.

In the years following the formulation of the dual porosity model for naturally fractured reservoirs, solutions of the coupled partial differential equations for the pore and fracture fluid pressures were obtained using finite difference techniques. While these simulations can provide accurate solutions, often in a complicated geometry covering the entire reservoir, the large number of computations involved made them cumbersome for analysis of large data sets. In response, an alternative, faster, method of solution was developed in the 1980s. For a simplified geometry, Laplace transform solutions were developed, in which the transformed solutions were inverted numerically, using, typically, the Stehfast algorithm.

Several authors have noted that analytic approximations can be developed for certain ranges of parameter values (referring to the Warren-Root dimensionless parameters defined below) appropriate for shale gas reservoirs. It will become apparent later in this paper that for typical shale gas reservoirs the interporosity flow coefficient (or transmissivity), $$\lambda$$ , is very small and this allows asymptotic approximations to be derived for the Laplace-transformed solutions. Since these models still require numerical inversion of the transformed solution, it would be more accurate to label them semi-analytic models. They have advantages over full numerical solutions in terms of speed of calculation and in the added value they bring to understanding the flow characteristics and the impact of the reservoir and completion parameters on production.

2.3. Development of new semi-analytic solution

Recently, we have taken the idea of developing asymptotic solutions one stage further. We have developed perturbation solutions for $$\lambda <<1$$ directly from the dual porosity partial differential equations, thereby removing the necessity for Laplace transforms altogether. The greater simplicity and enhanced understanding afforded by these solutions will become apparent as we proceed. (Full details are available in an internal EGI report [14].) The result is similar to the simple linear flow model that is currently gaining favor in the literature, but has some notable advantages:

• The model does not make the a priori assumption of linear flow into a sequence of transverse fractures.

• The model does not make the a priori assumption of infinite fracture conductivity, and allows an estimate to be made for the fracture pressure loss.

• Identification of the end of the linear infinite-acting flow period and development of the ensuing PSS solution is made explicit.

• Provision of a solution form that facilitates fast and easy production data analysis.

• Identification of the reservoir and completion parameters that are the greatest (primary) determinants of productivity.

• Solution scheme that permits rational extension to include other physical processes, such as desorption.

For simplicity we restrict attention in this paper to single-phase flow in the matrix and assume that gas is produced at constant bottom hole pressure; we shall also neglect the impact of gas desorption. For the purposes of the present discussion the most important part of the solution is the leading order solution for the reservoir pseudo-pressure, which satisfies a standard diffusion equation.

$$\frac{{\partial }^2{m}_{D{m}_0}}{{\partial }^{}{z}_D{}^2}=\frac{{\partial }^{}{m}_{D{m}_0}}{\partial t_D}$$ E1

Here we have used dimensionless parameters (denoted by subscript D); full definitions are provided in the Nomenclature section later in this paper. Dimensionless distance normal to the fracture face, $$z_D$$ , is scaled on L/2 (i.e. half the inter-fracture spacing), and dimensionless time, $$t_D$$ , is scaled on a time scale,

$$t_m=\frac{c{\phi }_m\mu {\left(\raisebox{1ex}{$L$}\!\left/ \!\raisebox{-1ex}{$2$}\right.\right)}^2}{k_m}$$ E2

that characterizes pressure diffusion in the matrix

The appropriate initial and boundary conditions are

$$\begin{array}{l}{m}_{D{m}_0}=0at{t}_D=0\\ m_{D{m}_0}=m_{D{f}_0}=1at{z}_D=0,\frac{\partial m_{D{m}_0}}{\partial z_D}=0at{z}_D=1\end{array}$$ E3

The leading order influx from the matrix into the fracture network is given by

$$q_{D0}={\left|\frac{{\partial }^{}{m}_{D{m}_0}}{{\partial }^{}{z}_D}\right|}_{z_D=0}$$ E4

At downhole conditions the mass flowrate from the matrix into the fracture network is

$$q_m=q_{ch}{q}_{D0}$$ E5

The dimensionless flowrate is defined in equation (4) and the characteristic mass flowrate is defined by

$$q_{ch}=A\frac{M{k}_m}{2R{T}_w}\frac{m_{ch}}{Z_{ch}}$$ E6

Here A denotes the productive fracture surface area. The total mass flow rate measured at surface conditions is

$$q_s=\frac{q_m}{{\rho }_s}=A\frac{k_m}{Z_{ch}}\frac{T_s}{T_w}\frac{Z_s}{p_s}{m}_{ch}$$ E7

Here subscript s denotes surface conditions. Note that $q_s$ is actually a volumetric flowrate and is measured typically in units such as scf/s.

The dimensionless flowrate defined in equation (4) may be readily calculated in terms of the dimensionless pseudo-pressure, either from the full numerical solution of the diffusion equation, (equation (1)), or from the early-time infinite-acting approximation to it. Both solutions provide very useful information and insights into the variation of the production rate with time. The dimensionless cumulative production is defined by

$$Q_{D0}={\displaystyle {\int }_0^{t_D}{q}_{D0}d{t}_D}$$ E8

The diffusion equation (1) is readily solved using standard numerical schemes as made available in mathematical software such as MATLAB. To complement this solution we have obtained an analytic approximation valid while the change in pressure has not been impacted by neighboring boundaries or fractures—often referred to as the infinite-acting approximation. The early-time approximation to the dimensionless inflow rate is

$$q_{D0}=\frac{1}{\sqrt{\pi t_D}}$$ E9

And the early-time approximation to the dimensionless cumulative production is

$$Q_{D0}=2\sqrt{\frac{t_D}{\pi }}$$ E10

Figure 5 compares the full and early-time solutions for the dimensionless flow rate and cumulative production against dimensionless time. For convenience we have used a log-log plot here. There are several major features of this plot that are worthy of further comment:

• The early-time (infinite-acting) solution provides a good approximation to the full numerical solution for dimensionless time, $$t_D$$ , less than about 0.15, which translates to about 3 years in dimensional terms.

• During this time frame, the flowrate is represented by a straight line of slope – ½ and the cumulative production is represented by a straight line of slope + ½.

• Boundaries (or in this case neighboring fractures) begin to influence the flow after this point in time and the solution departs from the simple linear dependence on the square-root of time. This is also evident in the numerical solutions presented by Bello and Wattenbarger [15] in their Figures 5, 7, 8 and 10 and in the field data shown in their Figure 1.

• The dimensionless cumulative production approaches a final value of 1 as it should because of the way we have defined the characteristic scales in our non-dimensionalization of the equations.

• The error incurred in estimating the cumulative production by extrapolating the infinite-acting solution beyond its region of validity is apparent from Figure 5.

• Almost 90% of the total gas that can be produced has been produced by the time $$t_D=1$$ . This then provides a simple interpretation of the matrix diffusion time as the time (in real terms) to produce 90% of the gas available.

Figure 6 shows the full and infinite-acting solutions for cumulative production plotted against the square-root of dimensionless time. As expected from equation (10) the early-time solution is well represented by a straight line with slope $\raisebox{1ex}{$2$}\!\left/ \!\raisebox{-1ex}{$\sqrt{\pi }$}\right.$ or 1.128. Later in this paper we develop this plot as the basis of our production data interpretation technique.

In anticipation of the application of these results to analysis of production data, it is useful to provide an expression for the cumulative production in dimensional terms. Analogous to equation (6) we define cumulative production at downhole conditions by

$$Q_m=Q_{ch}{Q}_{D0}$$ E11

The characteristic cumulative production scale is defined as

$$Q_{ch}=q_{ch}{t}_m=A\frac{M{k}_m}{2R{T}_w}\frac{m_{ch}}{z_{ch}}{t}_m$$ E12

The time scale $t_m$ is defined in equation (2) and $$Q_{D0}$$ is defined in equation (8).

Analogous to equation (9) we define cumulative production at surface conditions by

$$Q_s=\frac{Q_m}{{\rho }_s}=A\frac{k_m}{Z_{ch}}\frac{T_s}{T_w}\frac{Z_s}{p_s}{m}_{ch}{t_{m}Q}_{D0}=Q_{ch,\mathrm{ }s}{Q}_{D0}$$ E13

Figure 5.

A comparison of the full and early- time solutions for the dimensionless flow rate and cumulative production.

Figure 6.

A comparison of the full and early- time solutions for the dimensionless cumulative production.

In view of the wide applicability of the early-time solution, it is useful to state the form taken by equation (13) during the infinite-acting period. Using the early-time approximation given in equation (9), we see that

$$Q_s=A\frac{k_m}{Z_{ch}}\frac{T_s}{T_w}\frac{Z_s}{p_s}{m}_{ch}2{\sqrt{\frac{t_m}{\pi }}}_{}{\sqrt{t}}_{}$$ E14

If we now use the definition of $T_m$ (equation (2)), we can express the cumulative production at surface conditions as

$$Q_s=C_P\sqrt{t}$$ E15

Where

$$C_P=2A\frac{T_S}{T_W}\frac{Z_S}{p_S}{m}_{ch}\sqrt{\frac{c{\phi }_m{k}_m}{\pi \mu }}$$ E16

The coefficient $C_P$ represents the slope of the dimensional equivalent of the straight line in Figure 6 and is of fundamental importance in the subsequent development of this paper. For now we will observe that $C_P$ characterizes the early-time solution in a form that is easy to estimate from production data. We shall refer to it as the “Production Coefficient”. (We have adopted this terminology in recognition of the similarity of this result to a well-known expression for the leakoff rate of a compressible fluid from a fracture to a reservoir filled with the same fluid (see, for example, [15]).

3. Production data interpretation

The analysis described above suggests that for a substantial part of a shale gas well’s production history, the cumulative production varies linearly with the square-root of time. The coefficient $C_P$ represents the slope of the straight line in a plot of Q against $\sqrt{t}$ and characterizes the early-time solution in a form that is easy to compare with production data The time scale, $T_m$ , defined in equation (2) defines the upper limit of the applicability of the linear flow regime and allows us to characterize the production rate once boundary-effects have become important.

We illustrate this production analysis technique by comparing production from a group of wells in the Barnett shale. Figure 7 shows a conventional plot of production rate against time for several wells that had been producing for at least 5 years in an area of the Barnett field. The data was obtained from a public data base and we have plotted the production rates at yearly intervals. For clarity of presentation we have connected the data points by smooth lines. This “conventional” plot reveals nothing about the relative decline rates of the wells or provides much insight into the flow regime(s). The same data sets have been plotted in the new format in Figure 8. It is immediately apparent that for most of these wells the data falls on straight lines as expected from our analysis. The slope of these lines is readily measured and provides an estimate for the production coefficient, $C_P$ . Estimation of $C_P$ is quick and easy and provides us with a new metric with which we can quantify the productivity of these wells. Again, we explore these results in more detail later, but for now we note that the linear flow regime extends beyond at least 5 years, since there is no indication at this point of departure from the straight line in this plot. This fact alone sets some bounds on the fracture spacing and the matrix permeability.

In Figure 8 we have omitted the first year’s cumulative production data. Generally, early-time production data is quite severely impacted by variable drawdown conditions and so we should not expect a good straight line fit at that time. Analysis of this regime is discussed at length elsewhere [14].

Figure 7.

Production data from a group of Barnett wells plotted in the conventional format

Figure 8.

Same production data from Figure 7 plotted in the new format.

4. Identification of production drivers: Nature versus nurture

We have established the applicability of this new analysis technique to data from very many wells in many shale gas plays across the US and Canada, though there is neither space nor time available to discuss this in detail in the present paper. Following on from that analysis, we will now go on to discuss the results in more depth and begin to draw some tentative conclusions about the production drivers, by examining the parameters that together constitute the formula for $C_P$ in equation (15). We may divide these parameters broadly into two groups—those that characterize the nature of the reservoir and those that characterize how we “nurture” the reservoir. Specifically the parameters are:

Nature:

• Matrix permeability, $k_m$

• Matrix porosity, $${\phi }_m$$

• Gas viscosity, $$\mu$$

• Gas compressibility, c

• Initial reservoir pressure, $$p_i$$

• Reservoir temperature, T

Nurture:

• Bottom hole flowing pressure, $$p_w$$

• Productive fracture surface area, A

In developing these results we are constrained by the requirement that $\lambda$ <<1, where

$$\lambda =\frac{12k_m{r}_w^2}{L{c}_f}$$ E17

This requirement sets some bounds on the fracture network characteristics, but these are generally easily met for shale gas reservoirs. For given values of the matrix permeability and the wellbore radius, the combination of fracture spacing and fracture conductivity must be sufficiently large.

It is apparent that given these conditions production for a large part of the production history of these wells depends upon the parameters listed above. We note in particular that history matching production data over this flow period furnishes only one parameter and that is the production coefficient, $C_P$ . That is all. The square-root of time behavior is inherent to the physics of the flow: i.e. linear flow into a network of (effectively infinitely-conductive) fractures. It is not at all surprising that conventional history-matching techniques using reservoir simulators give non-unique answers: many different values of the parameters in the list above can together constitute the same value of the production coefficient. Moreover this formulation tells us what parameters have little effect on the history-matching process, including the precise value of the fracture conductivity. Even if history-matching were attempted in terms of dimensionless parameters, it is apparent that the result is insensitive to $\lambda$ provided that it is small enough.

We need to elaborate at this point on the parameter A defined above as the “productive fracture surface area”. This is the area of the fractures in contact with the reservoir that serve as the channels that convey gas from the matrix to the wellbore. We can make no assertion at this point about whether these fractures are natural fractures or propped or unpropped hydraulic fractures and nor can we say anything (yet) about their spacing or their lengths or indeed their number and location. All we can infer from the production data analysis is the total productive fracture surface area.

Some insight about the spacing of these productive fractures can be obtained by examining the time scale of pressure diffusion in the matrix. We demonstrated earlier that we may expect the root-time solution to be valid until neighboring fractures begin to compete with one another for production. In other words, until pressure diffusion in the matrix can no longer be considered to be independent of the fracture spacing. According to the analysis presented above we should expect the cumulative data to deviate from a straight-line in the root-time plot for t> 0.15 $$t_m$$ .

If we could detect the time at which this departure occurs then, we have some information with which to estimate the productive fracture spacing. Even if the entire production history to date is in the linear flow regime, we can make an estimate of a lower bound on the fracture spacing. As we see later, the fracture spacing is surprisingly large for typical shale gas plays.

We have now analyzed many shale gas production data sets using our proposed technique and have found the square-root fit to be very good. Based on this and on the mathematical analysis that supports that technique we have concluded that the production rate declines inversely with the square root of time. As we have discussed above, this is a consequence of the dominant production process of linear flow into a network of fractures. The decline rate is therefore fully determined by the physics of the production process. We should not expect to see any significant variation from well to well, from vertical to horizontal wells or indeed form play to play. What does vary is the multiplier, the production coefficient ${, C}_P$ , which as we have demonstrated elsewhere depends on many factors, principally the reservoir quality, the reservoir and bottom hole pressure and the productive fracture surface area.

4.1. Example: Barnett shale production data analysis

As we have indicated above, it is relatively straightforward to use this new technique to analyze production data whether it is on a well-by-well basis or averaged over a play or area within a play. In essence, the process consists of three steps:

1. From the daily (or monthly or yearly) production data calculate the cumulative production for each well at different points in time.

2. Plot cumulative production against the square-root of time.

3. Estimate the Production Coefficient form the slope of the best straight line fit to the data.

The Barnett shale is a good starting point for a more in-depth data analysis, since production data is readily available from public databases and, moreover, that data extend over many wells for long periods of time. The Barnett shale occupies several counties in North Texas. It is broadly bounded by geologic and structural features and may be divided according to estimates of maturity into a gas window and oil window. Historically the major development has been in a core area located to the north of Fort Worth (Figure 9), but more recently expansion has occurred to the south and to the west. To date many thousands of wells have been drilled and completed in the Barnett, initially vertical, but now almost entirely horizontal.

It has become common practice to sub-divide the Barnett play into three areas, described as the Core, Tier 1 and Tier 2. For convenience we may define these areas according to county as follows

• Core region: Denton, Tarrant and Wise counties, comprising 2974 horizontal wells and 3886 vertical wells

• Tier 1 region: Hood, Johnson and Parker counties, comprising 3865 horizontal wells and 251 vertical wells

• Tier 2 region: all other counties, comprising 687 horizontal wells and 401 vertical wells

It is apparent form this cursory division that the fraction of wells that were completed horizontally shifts from 43% in the Core area to 94% in the Tier 1 area, which reflects the development of technology with time and the spread of drilling with time to the outer areas. At the date of these figures (2009) Tier 2 was relatively unexploited.

The result of this detailed analysis (Figure 10) allows us to quantify the production variations in the Barnett Core, Tier 1 and Tier 2 areas and to distinguish the impact of horizontal and vertical well completions on the productivity. In a sense this represents a first, somewhat crude, pass at distinguishing the impact of nature (in the sense that reservoir properties depend on location, with the core area providing more fertile ground than Tier 1 or Tier 2) and nurture (in the anticipation that horizontal well technology provides more productive fracture surface area than does vertical well technology).

Figure 9.

Development of the Barnett shale in North Texas

In Figure 10 we have shown the cumulative distribution of production coefficient for each of the six categories defined above. The plots should be interpreted as follows. For each category the probability that a well has a specified value of the production coefficient in excess of the value on the x-axis can be read off the y-axis. For example the probability of a horizontal well in the Core having a production coefficient in excess of 0.75 (bcf/yr^0.5) is about 8%.

It is apparent that, as is to be expected, wells in the Core have better production characteristics than wells in Tier 1 and wells in Tier 2 and that in general horizontal wells have better production characteristics then vertical wells. It is interesting in this context to examine the variation of production coefficient in the core area in more detail. Figure 11 shows the location of ten of the wells in the core area with high values of the production coefficient (in green), 10 of the wells with medium values (in blue) and 10 wells with low values (in red). It is apparent that there appear to be sweet spots even within the core area, but there are substantial outliers and there are some relatively poor wells close to better wells.

Figure 10.

Distribution of values of the Production Coefficient for horizontal and vertical wells in the Core region, Tier 1 region and Tier 2 region of the Barnett shale.

Figure 11.

Preliminary identification of sweet spots in the core area of the Barnett. Each black dot represents a well; 10 wells with high Production Coefficients are identified in green, 10 medium performers are in blue and 10 poor performers are identified in red.

5. Implications and deductions

One of the advantages of the semi-analytic method outlined in this paper is that it enables us to make certain deductions about the magnitude of the parameters that drive the productivity of the well. In particular we can make some inferences about the magnitude of the fracture surface area through which the gas is produced and about the likely spacing of the productive fractures.

5.1. Productive fracture surface area

Our analytic solution allows us to relate the Productivity Coefficient $C_P$ to a group of parameters that may be roughly divided into those that characterize the nature of the reservoir and those that characterize the impact of the completion and stimulation strategy. In our formula for $C_P$ (equation (16)) perhaps the parameter that has the greatest uncertainty is the productive fracture surface area, A. This is the area of the fractures in contact with the reservoir that serve as the channels that convey gas from the matrix to the wellbore. We can make no assertion at this time point about whether these fractures are natural fractures or propped or unpropped hydraulic fractures and nor can we say anything (yet) about their spacing or their lengths or indeed their number and location. We can, however, make an estimate from the production data analysis of the total surface area of these productive fractures.

In Figure 12 we show estimates of the productive fracture surface area for typical strong, medium and weak performing wells in the core area of the Barnett in terms of values of the matrix permeability. In this calculation we have made reasonable estimate of the other parameters that impact the productivity coefficient such as the gas viscosity and compressibility and the matrix porosity.

Several features of this plot merit discussion:

• Well productivity increases with productive fracture surface area and with matrix permeability, as expected.

• Wells in this group typically have matrix permeabilities in the range 100-300 nd and this implies that the productive fracture surface area is in the range 1-6 million square feet (Msqft): the higher the permeability, the less fracture surface area is needed to achieve the given productivity.

• If the matrix permeability was as low as even 10 nd, the required fracture surface area would approach 100 Msqft. On the other hand, matrix permeabilities of the order of 1 μd would require less than about 1 Msqft of productive fracture surface area.

Figure 12.

Estimate of productive fracture surface area for specified matrix permeability. The three curves were developed based on analysis of Barnett shale well data using the EGI semi-analytic production model.

To place these numbers in context we note that a fracture of height 200 ft and half length 200 ft has surface area of 0.16 Msqft. Thus, 20 of these fractures would have a fracture surface area of 3.2 Msqft, which is a perfectly plausible estimate of the hydraulic fracture surface area created with modern multi-stage fracturing techniques. (Note that 20 such fractures would be spaced about 150 ft apart in a 3000 ft lateral.)

Figure 12 was developed on the basis of production data from wells in the core area of the Barnett shale, but the results apply, at least qualitatively, to other shale or tight gas plays. For example for more conventional tight gas plays for which the permeability is of the order of 1 μd or more, we should expect respectable productivity with only one such hydraulic fracture, which reinforces our experience that a vertical well with a single bi-wing fracture may be adequate for those reservoirs, but not for shale gas plays. Conversely the productive fracture surface area for economic production from ultra-tight shale plays (such as the shallow shale plays described earlier in this paper) cannot be achieved by producing from the hydraulic fractures alone.

We have been careful so far to make no formal distinction between the natural fractures and the hydraulic fractures in so far as productivity is concerned. All we have demanded is that their conductivity is sufficiently large that $$\lambda <<1$$ , which should not in principle present too great a restriction on the fracture conductivity whether it is associated with propped fractures or unpropped fractures. On the basis of our data analysis we should expect a productive fracture surface area in the range of 1-6 Msqft. How then is that area created and what are the implications of this figure?

A crude estimate of the fracture surface area that is created by pumping large volumes of frac fluid may be made by performing a mass balance and assuming that none of the fluid has leaked off or imbibed into the formation over the time in which the fracture network is created. For a total fluid volume V and assuming a created average fracture width w during pumping, the total fracture surface area may be estimated at

$A=2\frac{V}{w}$

(using any consistent set of units of course). If, for example, 100 million gallons of frac fluid were pumped and the assumed frac width was 0.2 inches, then the total surface area would be about 100 Msq ft. Clearly, this is far in excess of our estimate of the productive fracture surface area and would suggest that less than 10% of the created fracture surface area is actually productive. Naturally, this raises all sorts of other questions concerning the efficiency of this process, which we plan to address in a future project.

A similar mass balance for the proppant placed in a typical job enables an estimate to be made of the surface area of propped fractures. If we make some estimate of the likely width (0.1 in) and porosity (0.4) of a propped fracture (after closure), it appears that a propped fracture surface area of the order of a few million square feet is quite plausible.

5.2. Productive fracture spacing

Some insight about the spacing of these productive fractures can be obtained by examining the time scale of pressure diffusion in the matrix. We demonstrated earlier that we may expect the root-time solution to be valid until neighboring fractures begin to compete with one another for production. In other words until, pressure diffusion into a fracture can no longer be considered to be independent of the fracture spacing. According to the analysis presented earlier we should expect the cumulative production data to deviate from a straight-line in the root-time plot for $$t>0.15T_m$$ . If we could detect the time at which this departure occurs then, we have some information with which to estimate the productive fracture spacing. Even if the entire production history to date is in the linear flow regime, we can at least make an estimate of a lower bound on the fracture spacing.

It is instructive to estimate the matrix diffusion time for typical values of the fracture spacing and matrix permeability. The results are shown in Figure 13. The diffusion time increases quadratically with the fracture spacing and inversely with the matrix permeability. Typical values for the diffusion time are quite low. For example, if, as we expect, linear flow continues for at least 3 years, then we should expect to see a diffusion time of the order of 20 years. Figure 13 suggests that the productive fracture spacing is likely to be of the order of 100 ft or more.

Figure 13.

The impact of fracture spacing on the time to produce 90% of the gas in place.

In figure 13 we have identified the diffusive time scale with the matrix drainage time. As we showed above, 90% of the total gas in the pore space between the fractures has been drained by this time. A time scale of about 20 years is at least consistent with the industry estimates of the effective production lifetime of these wells. It is worth noting here the consequences of much smaller fracture spacing. For a fracture spacing of only 10 ft, we estimate that 90% of the total gas production will have occurred within the first few months of production, which is quite unrealistic. Note also that the surface area of planar fractures only 10 ft apart in a 3000 ft lateral would be of the order of 150 Msq ft, which again is unreasonably large.

6. Conclusions

A common view of production mechanisms in shales is “because the formations are so tight gas can be produced only when extensive networks of natural fractures exist” [6]. To this extent gas production from some of the shallower (Devonian) shales is similar to gas production from coal. As we have discussed earlier in this paper, we expect that the deeper gas shales differ in this respect.

Using a new semi-analytic production model, we have analyzed production data from a number of shale gas wells in several different North American shale gas plays. Interpretation of the results suggest that productivity is largely determined by a small group of parameters that may be decomposed into two sub-groups representing the nature of the reservoir (such as matrix permeability and porosity) and what we may term our (engineering) attempts at nurture (including completion and stimulation parameters). Of key importance is the productive fracture surface, which unfortunately is difficult to estimate a priori. However, our interpretation of the production data suggest the following

• Productive fracture surface area ~1-6 Msqft and probably within 2-4 Msq ft.

• The volume of these productive fractures is very much less than the volume of water pumped, but

• Productive fracture volume scales approximately with the volume of proppant placed.

• Typically, there is no indication of fracture interference during production even after several years, which suggests that the productive fracture spacing is at least 100 ft.

• Time to drain 90% of the fractured region or matrix blocks: ~10-20 years

We are led to the conclusion that almost all the fracturing fluid pumped during a multi-stage horizontal well fracturing operation in the shales serves to open a vast, and possibly complex, network of natural fractures and that these fractures do not make a significant contribution to the well’s productivity. We are led inevitably to questions concerning the conductivity of these, largely unpropped, fractures and to investigate the rock and fluid mechanisms that seemingly prevent them from being productive. The role of the fracturing fluid (usually slickwater) in this process should now be investigated from this new perspectivel

Nomenclature

k – permeability

$${\phi }_{}$$ – porosity $$\mu$$ – gas viscosity $\lambda$ – dual porosity transmissivity factor (defined in equation (17))

c – gas compressibility

cf – fracture conductivity

$$p$$ – pressure

m – gas pseudo-pressure

$C_P$ – production coefficient (defined in equation (16))

q – production or flow rate

Q – cumulative production

r – radius (wellbore)

t – time

T – reservoir temperature

A – productive fracture surface area

L – fracture spacing

z – co-ordinate normal to fracture surface

Z – real gas compressibility factor

Subscripts

ch – characteristic

D – dimensionless

m – matrix

i – initial

s – surface

w – wellbore