International Journal of Clean Coal and Energy
Vol.04 No.02(2015), Article ID:54685,10 pages

Modeling and Simulation of Natural Gas Production from Unconventional Shale Reservoirs

Gary Feast1, S. Sina Hosseini Boosari2, Kim Wu1, John Walton3, Zufang Cheng3, Bao Chen3

1Novel Energy, Inc., Pittsburgh, PA, USA

2Petroleum and Natural Gas Department, West Virginia University, USA

3Subsurface Shale Group, Columbus, OH, USA


Copyright © 2015 by authors and Scientific Research Publishing Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY).

Received 22 February 2015; accepted 12 March 2015; published 17 March 2015


Modeling and simulation of unconventional reservoirs are much more complicated than the conventional reservoir modeling, because of their complex flow characteristics. Mechanisms, which control the flow in the reservoir, are still under the investigation of researchers. However, it is important to investigate applications of mechanisms which are present to our knowledge. This paper presents the theory and applications of flow mechanisms in unconventional reservoir modeling. It is a well-known fact that most of the reservoir flow problems are non-linear due to pressure dependency of particular parameters. It is also widely accepted that fully numerical solutions are costly both computational and time wise. Therefore, the presented model in this paper follows semi-analytical solution methods. Gas adsorption in unconventional reservoirs is the major pressure dependent mechanism; in addition existence of natural fractures is also taken considerable attention. This paper aims to investigate combined effect of existence of natural fractures gas adsorption, and gas slippage effect while keeping the computational effort in acceptable range. Unlike the existing literature (Langmuir is widely used), BET multi-layer isotherm employed in this paper for gas adsorption modeling. A modified dual porosity modeling is used for natural fracture and gas slippage effect modeling. For model verification purposes a history matched is performed with real field data from Marcellus shale. The proposed model in this paper shows a good agreement with the field data. It is observed that BET isotherm models early time production performance more accurately than Langmuir isotherm. It is also concluded that gas adsorption significantly improves the production performances of unconventional reservoirs, with natural fractures. In addition, gas slippage has a slight effect in long term production.


BET Desorption Isotherm, Shale Gas Reservoirs, Hydraulic Fracturing

1. Introduction

Technological improvements in the oil and gas industry made unconventional reservoirs a hot hydrocarbon resource especially in the last decade. Different from conventional resources, unconventional reservoirs such as shale rocks have some distinct characteristics. First and most importantly shale reservoirs are the source rock with their organic content. It means that unconventional reservoirs are both source and cap rock for hydrocarbons. Secondly, unconventional reservoirs are extremely tight with nano-Darcy scale permeabilities. In addition, this extremely tight condition also alters the conventional pore sizes. Pore size scale is an important fact, because it is commonly known that Darcy-flow is the dominant flow regime in conventional reservoirs. However, when it comes to shale resources, other flow mechanisms such as slippage flow might be important and needed to be taken into consideration, due to small scale of pore size distribution. Thirdly, not only hydraulic fractures, but also natural fractures are the key elements of unconventional reservoirs which might play a crucial role in production. This paper aims to address, these flow mechanisms in a rigorous mathematical flow model, and highlights their importance. In addition, a history match was performed with newly developed reservoir model for verification purposes.

As it is stated earlier, controlling flow mechanisms in unconventional reservoirs is quite complex and needed to be considered in reservoir flow modeling. For this purpose, Langmuir isotherm is a preferred isotherm model to model gas desorption. Recently, researchers [1] claimed that Brunauer, Emmet and Teller (BET) isotherm can predict some of the shale reservoirs desorption profile better than Langmuir isotherm. Moving from gas adsorption, gas slippage effect is also discussed by researchers and various mathematical models are presented to our knowledge. Scientist [2] modified permeability approach and [3] [4] apparent gas slippage term was commonly used in the literature for gas slippage effects. Natural fractures are another phenomenon both for conventional and unconventional reservoirs. Dual porosity models in analytical approaches are customarily used for natural fracture modeling. Although, fundamentals of dual porosity models are coming from conventional reservoirs, many dual porosity models specific to unconventional reservoirs are proposed by researchers in the literature [5] - [10] . Until recently [11] , researchers mainly focused on one of these control mechanisms rather than combining them in a compact model. In this paper, we proposed our semi-analytical reservoir model with visiting the theory and applications of different flow mechanisms.

2. Derivation of New Flow Model for Unconventional Reservoirs

Apart from commercial reservoir simulators, analytical models include fundamental physical solutions, and practical tools for reservoir modeling. However, when the non-linearity comes into the picture, analytical models are not adequate enough. Therefore, it is necessary to implement some degree of numerical methods into reservoir flow solution. In order to minimize the necessary input parameters and computational algorithm, our solution follows semi-analytical approach. Having said that, the only steps which are involving numerical methods are gas desorption and gas slippage terms in our formulation. These two terms create non-linearity, because of their pressure dependency. In order to overcome from these issues, a methodology proposed by [12] - [14] is followed in our solutions. The average reservoir pressure calculation and using the calculated average reservoir pressure for the pressure dependent variables calculation is the main algorithm of the approach. Gas diffusivity equations used in our model are stated starting from Equation (1), associated mechanisms which are considered in our flow model such as dual porosity model, gas adsorption model, and slippage flow are given in further details.

( ρ g u g ) x + F ( L , t ) = [ ϕ ρ g + ( 1 ϕ ) ρ g s c G ρ b ] t (1)

where ug represents the velocity of gas, ρg shows free gas density; ρgsc is the gas density in standard condition, G represents potential releasable-gas content in scf/ton in particular model, ρ b is shale matrix density and ϕ denotes rock porosity. F represents the source term which represents the mass influx from matrix to natural fractures per unit time step. The details of the natural fracture flow modeling are given under the related section. L is the characteristic length, here in this study it is considered as hydraulic fracture half-length. Details of the gas adsorption model and slip flow phenomena are given under the related sections. It must be noted that not only the general form of flow equation, but also different mechanisms include pressure dependent terms and increase the order of non-linearity. Equation (1) is highly non-linear because of the pressure dependent parameters, however combining pressure dependent parameters under the same expression and calculating the average reservoir pressure for each time step, and using the average pressure as a reference point to update the pressure dependent parameters, we were able to manage computational effort and robustness of the our reservoir model.

Our new model considers evenly distributed hydraulic fractures along the horizontal well. The solution obtained from our solutions considers one representative hydraulic fracture and its corresponding natural fractures. Gas adsorption in our formulation can be applicable for different methods, here in this paper we used Langmuir and BET isotherms. Slippage flow is also taken into account for extremely small pore size nature of shale reservoirs.

3. Gas Adsorption in Unconventional Reservoirs

Coming from the nature of shale gas reservoirs, there are mainly three sources of natural gas. These are free gas, adsorbed gas on the organic material, and the organic material (kerogen most of the case) itself. Once the production starts from shale resources, first the available free gas will be produced. Therefore, there will be no equilibrium in the pore space and the adsorbed gas starts desorbing. Eventually, new gas molecules from the organic content inner pores start migrating to the surface of the kerogen and releases from there as free gas [15] .

Langmuir isotherm is customarily preferred to model gas adsorption process and has been placed in the fundamental gas diffusivity equation as a gas source term which is pressure dependent (creates non-linearity in the diffusivity equation). Equation (2) represents Langmuir equation.

G = V L P P + P L (2)

where G is the potential releasable gas amount, P represents the average reservoir pressure, VL stands for Langmuir volume constant and PL is Langmuir pressure constant. These constant are formation specific, and can be determined by laboratory experiments.

In addition to Langmuir adsorption model, there are other models available to our knowledge. However, until recently laboratory experiments were limited and their applicability needed to be approved scientifically by the researchers. For this purpose BET [16] isotherm is studied by [17] .

Differently from Langmuir isotherm BET assumes multilayer adsorption phenomena instead of single layer. Hence, the model predicts more adsorb gas production from unconventional organic rich reservoirs. Starting from Equation (3) governing equations for BET model are given below.

G = v m C p ( p 0 p ) [ 1 + ( C 1 ) p / p 0 ] (3)

where p0 represents the original gas pressure, vm denotes maximum adsorb gas volume for a single layer, and C symbolizes a heat constant which is defined in Equation (4).

C = exp ( E 1 E L R T ) (4)

where E1 and EL are the first layer and higher degree layers of adsorption heat. R is the constant, and T represents the reservoir temperature.

In this paper, we used both of Langmuir and BET isotherm models with our semi-analytical model in order to compare their ability to predict unconventional reservoir production performance.

4. Slippage Flow and Pressure Dependent Natural Fracture Conductivity in Unconventional Reservoirs

Since the pore sizes are extremely low in shale reservoir environment it is important incorporate slippage flow in our theoretical reservoir modeling. Nanometer order of reservoir matrix pore sizes dictates the controlling flow mechanism as slippage flow which is pressure dependent. To mathematically express this phenomena, apparent permeability approach is proposed. Slippage flow is commonly referred in the literature as Klinkenberg effect. For this purpose Ertekin et al. (1986) proposed an apparent gas slippage term which is given in Equation (5).

b a = D g μ g c g p k (5)


D g = 31.57 M g k 0.67 (6)

µ is the gas viscosity, c is the compressibility of gas, P reservoir pressure at current time step, k represents the permeability, and Mg is molecular weight of the gas. Finally, gas slippage term can be replaced in the general form of Klinkenberg (1941) function which is given in Equation (7).

k a p p = ( 1 + b a p ) (7)

where kapp represents the pressure dependent apparent permeability. As can be clearly seen in apparent permeability calculation is function of pressure, therefore creates a non-linearity in the flow equation. Slippage flow considered for the flow from matrix to natural fractures.

5. Natural Fracture Modeling for Unconventional Reservoirs

Unconventional shale reservoirs are hydraulically fractured, and this process opens the cemented natural fractures in the tight formations. Hydraulic fractures are the main conductive paths for reservoir flow. In addition, secondary flow channels which have relatively low conductivity then the hydraulic fractures are became open with hydraulic fracturing. However, the density of natural fractures is considerably high compare to hydraulic fractures, and needed to be considered in the reservoir model.

Dual porosity models are considered in unconventional reservoir modeling to model natural fractures. For instance [18] used classical dual porosity models [19] in their analytical model. Later, [7] proposed their dual porosity model which is specific to unconventional reservoirs, and takes Knudsen flow and natural fracture closure effects into account. They have developed their model with considering spherical matrix blocks. However, we derived same equations with considering slab matrix in linear coordinates. Our dual porosity scheme also considers slippage flow from matrix to natural fractures. Governing equations for the dual porosity modeling is given starting with Equation (8). Here the aim is to derive a transfer function which models flow from matrix to natural fractures. As it is a well know phenomena, pressure dependency of gas creates non-linearity, therefore we used pseudo gas pressure definition by following [10] . In addition, Laplace theorem is used to solve partial differential mathematical expressions with following [20] technique, and [21] numerical conversion algorithm is used to obtain current time step production rate.

In order to mathematically express the flux from the matrix to natural fractures defines the fundamental version of slab matrix transfer function. In our model we modified their equation by following [7] approach with implementing all the pressure dependent terms under the pseudo-pressure term and defining a pressure dependent term β to represent flow from matrix to natural fractures:

f ( L m D , s ) = s β [ 1 + ω ˜ λ ˜ 3 s tanh ( 3 ω ˜ s λ ˜ ) ] (8)

where f function represents the flow from matrix block to fractures, ω represents the transmissivity, λ represents the transmissivity of the matrix blocks and corresponding equations are given in Equations (9) and (10).

λ ˜ = k m L 2 k f i h f (9)

ω ˜ = ϕ m c t m 2 ϕ f c f (10)

In order to represent slip page flow in matrix blocks modified pseudo-pressure is given in Equation (11).

m m ( p ) = 2 p * P ( 1 + b a m p m ) p d p μ z (11)

For natural fracture pseudo-pressure:

m f ( p ) = 2 p * P p d p μ z (12)

To overcome the pseudo-pressure inequality between matrix and natural fractures we followed Ozkan et al. (2010) and Aybar et al. (2014) and defined the term β and included that term in the transfer function. This term is also pressure dependent and we used average reservoir pressure for the current time step to iteratively obtain the β (Equation 13).

β = m m ( p ) m f ( p ) (13)

6. Model Verification

It is necessary to verify our proposed model with an actual field data. As it is given under the gas adsorption section recent studies showed that Marcellus Shale samples can be modeled by BET adsorption model more accurately than the Langmiur adsorption model. Therefore, we used the available field data from Marcellus shale in order to verify our semi-analytical model. We also used the available experimental data for gas adsorption and its fitting parameters from the recent study by Eshkalak et al. (2014), and incorporated with our reservoir model. It is assumed that reservoir is naturally fractured. It is also assumed that the well is evenly hydraulically fractured, and producing from a single shale layer. Also, hydraulic fracture height and the producing shale thickness assumed to be equal. Three years of actual production data was available, and we have used that data to verify our semi-analytical model.

Table 1 gives the reservoir model details used in our simulations, and Table 2 shows the details of Langmuir

Table 1. Reservoir parameters for history matching.

Table 2. BET and Langmuir isotherms fitting parameters used in this study.

and BET adsorption models fitting parameters. We used the gas flow rate as a matching parameter with keeping bottom hole pressure as constant. For all the pressure dependent terms in our reservoir model, we calculated the average reservoir pressure for each time step. Consequently, we used these average pressures to update our pressure dependent terms for the corresponding time steps. This methodology first proposed by Aybar et al. (2014), we have extended their methodology by incorporating different type of adsorption parameters.

Figure 1 shows our history matching with the different conditions. As it can be clearly seen, taken gas desorption into account significantly increased the production rates. In addition, BET isotherm gives better history match with the field data compared to Langmuir isotherm. This results in agreement with the results of the current studies about Marcellus shale desorption modeling. Upon completion of history matching we have also verified our new semi-analytical model and continued on the parameter analysis in the next section.

7. Parameter Analysis

The aim of this study is proposed a new reservoir model with incorporating all the flow mechanisms. We verified our model with actual field data. Moreover, it is also essential to identify different mechanisms effect on production in long-term. For this purpose the history matched model is employed in our parameter sensitivity analysis, since it perfectly represents in-situ conditions. Langmuir and BET adsorption models, gas slippage and existence of natural fractures effects on cumulative production are evaluated independently in this section. 25 years of production time is selected for each condition, and our newly developed semi-analytical model is used to predict cumulative production performance. Table 3 includes the parameters and their ranges used in the parameter analysis.

The results for sensitivity analysis for BET and Langmuir isotherm are given in Figure 2. It is concluded that both BET and Langmuir gas adsorption considerations creates around 13% more gas production in 25 years. Nevertheless, as it is concluded in the model verification section BET model predicts Marcellus shale production behavior more accurately than Langmuir isotherm especially for the early production period. However this effect diminishes for the long term, because eventually the cumulative production will be controlled by the gas in place rather than the flowing mechanisms such as adsorption.

It is concluded that slippage flow effects on cumulative production performance is around 1%, which can be ignored in order to reduce the non-linearity of the governing flow equation. The results for the case considered in this study are given in Figure 3 for 25 years of production. This conclusion can provide some extended of understanding about slippage flow effects on production in longer term production scenarios.

The existences of natural fractures existence in unconventional shale reservoirs are extremely important. For the case studied in this parameter studies section which has 20 natural fractures with 0.05 md-ft conductivity, the cumulative production differences between homogenous (without natural fractures) and naturally fractured case is around 20% for 25 years of production shown in Figure 4. It is observed that the homogenous case cannot reach the absolute production limits due to extremely low matrix permeability condition. However, natural fractures enhance the effective permeability of the system and resulted in higher production performances for the same production time.

8. Summary and Conclusions

We presented our new semi-analytical reservoir model with incorporating fundamental flow mechanisms ob-

Figure 1. History match with the actual Marcellus shale field data.

Figure 2. BET and Langmuir isotherm parameter analysis results.

Figure 3. Cumulative gas production with and without slip- page flow.

served in shale reservoirs. Flow mechanisms considered in this study (gas desorption, flow in natural fractures, and slippage flow) are pressure dependent and create non-linearity in flow equations. Therefore, we have calculated average reservoir pressure numerically, and updated our pressure dependent parameters and obtained our simulation results accordingly. In addition, two different gas adsorption models are considered and studied in this study, and differences between these two models are underlined.

Conclusions from this study can be categorized as following:

Figure 4. Cumulative gas production with and without natural fractures.

Table 3. Data used for parameter analysis.

1) The presented model can accurately model the flow behavior of the shale resources.

2) A good history match is obtained for the Marcellus shale field data using the new semi-analytical model.

3) BET adsorption isotherm can be applicable for Marcellus shale especially for early production period.

4) Both BET and Langmuir isotherms significantly contribute to the production performance of the case under the study in this paper.

5) Existence of natural fractures resulted in 20% higher cumulative production within 25 years of production.

6) Slippage flow effect on production performance is extremely low, and can be ignored to reduce the non-li- nearity of the governing flow equation.

7) When all the flow mechanisms are considered together, it is observed that gas adsorption and existence of natural fractures are extremely important and dominate the flow in the reservoir to hydraulic fractures.

8) All fluid flow simulations are time-consuming specially in unconventional with unique features, hence this progress has intensive application in computational fluid dynamic using AI.

In addition to commercial reservoir models, this paper aimed to show that fundamental physics and mathematical solutions can accurately model the flow in shale reservoir with incorporating different flow mechanisms and has shown promising results in modeling the hydrocracking process. In this regard, the approach covers the full span of E&P industry.


  1. 1. Eshkalak, M.O., et al. (2014) Simulation Study on the CO2-Driven Enhanced Gas Recovery with Sequestration versus the Re-Fracturing Treatment of Horizontal Wells in the U.S. Unconventional Shale Reservoirs. Journal of Natural Gas Science and Engineering, 21, 1015-1024.

  2. 2. Ozkan, E., Raghavan, R.S. and Apaydin, O.G. (2010) Modeling of Fluid Transfer from Shale Matrix to Fracture Network. SPE Annual Technical Conference and Exhibition, Florence, 19-22 September 2010, Paper SPE 134830.

  3. 3. Aybar, U., Eshkalak, M.O., Sepehrnoori, K. and Patzek, T.W. (2014) Long Term Effect of Natural Fractures Closure on Gas Production from Unconventional Reservoirs. SPE Eastern Regional Meeting, Charleston, 21-23 October 2014, Society of Petroleum Engineers, Paper SPE 171010.

  4. 4. Eshkalak, M.O., Aybar, U. and Sepehrnoori, K. (2014) An Integrated Reservoir Model for Unconventional Resources, Coupling Pressure Dependent Phenomena. Eastern Regional Meeting, Charleston, 21-23 October 2014, Paper SPE 171008.

  5. 5. Perdomo, J.F.A., Civan, F., Devegowda, D. and Sigal, R.F. (2010) Accurate Simulation of Shale Gas Reservoir. Annual Technical Conference and Ehibition, Florence, 19-22 September 2010, Paper SPE 135564.

  6. 6. Eshkalak, M.O., Mohaghegh, S.D. and Esmaili, S. (2013) Synthetic, Geomechanical Logs for Marcellus Shale. Digital Energy Conference and Exhibition, The Woodlands, 5-7 March 2013, Paper SPE 163690

  7. 7. Ozkan, E., Brown, M.L., Raghavan, R. and Kazemi, H. (2011) Comparison of Fractured-Horizontal-Well Performance in Tight Sand and Shale Reservoirs. SPE Reservoir Evaluation & Engineering, 14, 248-259.

  8. 8. Eshkalak, M.O., Aybar, U. and Sepehrnoori, K. (2014) An Economic Evaluation on the Re-Fracturing Treatment of the US Shale Gas Resources. Eastern Regional Meeting, Charleston, 21-23 October 2014, Paper SPE 171009.

  9. 9. Boosari, S.S.H, Aybar, U. and Eshkalak, M.O. (2015) Carbon Dioxide Storage and Sequestration in Unconventional Shale Reservoirs. Journal of Geoscience and Environment Protection, 3, 7-15.

  10. 10. Aybar, U., Yu, W., Eshkalak, M., Sepehrnoori, K. and Patzek, T.W. (2015) Evaluation of Production Losses from Unconventional Shale Reservoirs. Journal of Natural Gas Science and Engineering, 23, 509-516.

  11. 11. Mengal, S.A. and Wattenbarger, R.A. (2011) Accounting for Adsorbed Gas in Shale Gas Reservoirs. The SPE Middle East Oil and Gas Show and Conference, Manama, 25-28 September 2011, Paper SPE 141085.

  12. 12. Eshkalak, M.O., Al-Shalabi, E.W., Aybar, U. and Sepehrnoori, K. (2014) Enhanced Gas Recovery by CO2 Sequestration versus Re-Fracturing Treatment in Unconventional Shale Gas Reservoirs. Abu Dhabi International Petroleum and Exhibition and Conference, in Abu Dhabi, 10-13 November 2014, Paper SPE 172083.

  13. 13. Langmuir, I. (1918) The Adsorption of Gases on Plane Surfaces of Glass, Mica and Platinum. Journal of the American Chemical Society, 40, 1403-1461.

  14. 14. Boosari, S.S.H., Eshkalak, M.O. andAybar, U. (2015) The Effect of Desorption-Induced Porosity-Permeability Changes and Geomechanics on Production from U.S. Shale Gas Formations. 49th U.S. Rock Mech. Symp.

  15. 15. Tao, Q., Ghassemi, A. and Ehlig-Economides, C.A. (2010) Pressure Transient Behavior for Stress-Dependent Fracture Permeability in Naturally Fractured Reservoirs. CPS/SPE International Oil and Gas Conference and Exhibition, Beijing, 8-10 June 2010, SPE 131666.

  16. 16. Aybar, U. (2014) Investigation of Analytical Models Incorporating Geomechanical Effects on Production Performance of Hydraulically and Naturally Fractured Unconventional Reservoirs. MSc Thesis, The University of Texas at Austin, Austin.

  17. 17. Patzek, T.W., Male, F. and Marder, M. (2013) Gas Production in the Barnett Shale Obeys a Simple Scaling Theory. Proceedings of the National Academy of Sciences of the United States of America, 110, 19731-19736.

  18. 18. Omidvar Eshkalak, M. (2013) Synthetic Geomechanical Logs and Distributions for Marcellus Shale. MSc. Thesis, West Virginia University, Morgantown.

  19. 19. Cho, Y., Ozkan, E. and Apaydin, O.G. (2013) Pressure-Dependent Natural-Fracture Permeability in Shale and Its Effect on Shale-Gas Well Production. SPE Reservoir Evaluation & Engineering, 16, 216-228.

  20. 20. Aybar, U., et al. (2014) The Effect of Natural Fracture’s Closure on Long-Term Gas Production from Unconventional Resources. Journal of Natural Gas Science and Engineering, 21, 1205-1213.

  21. 21. Eshkalak, M.O., Mohaghegh, S.D. and Esmaili, S. (2014) Geomechanical Properties of Unconventional Shale Reservoirs. Journal of Petroleum Engineering, 2014, Article ID: 961641.


BHP Bottom hole pressure, psi

ɷ Storativity

s Laplace transform parameter

m(p) Pseudo pressure, psi2/cp

mmD Dimensionless matrix pseudo pressure

mfD Dimensionless natural fracture pseudo pressure

Vm Matrix volume, ft3

Vf Natural fracture volume, ft3

P Pressure, psi

C Concentration constant related to the adsorption heat

Mg Molecular weight of the gas

PL Langmuir pressure, psi

VL Langmiur volume, ft3

km Matrix permeability, md

kf Natural fracture permeability, md

L Characteristic length, ft

hf Natural fracture thickness, ft

ʎ Transmisivity

ɸm Matrix porosity

ɸf Natural fracture porosity

df Rock characteristic parameter, psi−1

ctm Matrix total compressibility, psi−1

ctf Natural fracture total compressibility, psi−1