You are here
Home > BLOG > Flow Mechanics > Modeling and simulation of gas flow behavior in shale reservoirs

Modeling and simulation of gas flow behavior in shale reservoirs

Fig. 21 3-D dynamic pressure depletion.


Shale is a growing prospect in this world with decreasing conventional sources of fossil fuel. With the growth in demand for natural gas, there is impending need for the development of the robust model for the flow of shale gas (Behar and Vandenbroucke in Org Geochem, 11:15–24, 1987). So the major driving force behind the working on this major project is the unavailability of desired models that could lead to enhanced production of these wells and that too efficiently. This model mainly includes the movement of shale gas from tight reservoir through the conductive fractures to wellbore and production model of the decline in pressure inside the reservoir with respect to time.


Vamsi Krishna Kudapa, Pushpa Sharma, Vibhor Kunal, D. K. Gupta

Petroleum Engineering Department, University of Petroleum and Energy Studies, Dehradun, India

Received: 6 September 2016 / Accepted: 5 January 2017 / Published online: 2 February 2017
© The Author(s) 2017

This result has been further compared with the help of MATLAB so as to obtain a complete pressure-derived model. The result shows the applicability of this in the real-life projects where it is difficult to model the fractures and obtain the flow rate with them in fractures and how to set the production facilities becomes a question.


Shale is known as fine-grained, clastic sedimentary rock. The molecule size of shale is little which makes the interstitial spaces likewise little. Indeed, they are minute to the point that oil, regular gas and water experience issues traveling through the development. Shale can hence serve as a compelling top rock for oil and common gas (Firoozabadi 2012). Despite the fact that the interstitial spaces in shale are minute, they can take up a huge volume of the arrangement rock. This lets the shale to hold noteworthy measures of water, gas or oil and not have the capacity to adequately transmit them as a result of its low permeability. The petroleum business has beat these confinements of shale developments by utilizing level penetrating and hydraulic cracking to make build porosity and permeability inside the stone (Bustin et al. 2008).

Shale gas will be gas that is actually present in shale rocks. Sandstone rocks are known for high permeability, and gas can stream effortlessly through the stone. Interestingly, shale shakes for the most part have low permeability (Bustin et al. 2008).

Shale gas is viewed as an alleged “unusual gas,” together with “tight gas” with low permeability and “coal-bed methane” (CBM). While both traditional and capricious stores contain normal gas, it is the more intricate generation strategies that recognize the ordinary and offbeat store (Gong et al. 2011). Hydraulic breaking is regularly connected to capricious normal gas stores. India has immense stores of shale gas. As indicated by the accessible sources, India has around 300–2100 tcf evaluated gas setup in Indian shale gas bowls which is much bigger than stores that are accessible in Krishna–Godavari (D 6) Basin (Swami et al. 2013).

This paper mainly discusses about the modeling of gas flow from the matrix to the wellbore. The representation of the reservoir model includes a cube as a porous media, i.e., it contains pore spaces in which free gas is stored and also the adsorbed gas. Now, the gas in the cube (both free gas and adsorbed gas) will start flowing out inside the matrix to the fractures (induced). Many of these cube representations are put together and connected to the well bore.

In this paper, we have considered a updated dual-mechanism model. One porosity is the combination of matrix and natural fracture, and the second porosity is the hydraulic fracture. For this model, a nonlinear PDE equation has been developed which is then compiled using MATLAB to develop a simulator for calculating the shale gas production, by considering the matrix as a source term. The production data that are obtained from this model will describe the unique characteristics shale gas reservoirs.

A three-dimensional shale gas reservoir model was created. Three flow mechanisms (Darcy flow and non-Darcy flow) as well as gas adsorption and desorption mechanism were considered in this model. The flow in the matrix is considered as single-phase flow, and the production from this reservoir model is estimated for a period of 3 years and the results are validated by CMG-IMEX software.

Back ground literature

Modeling of unconventional gas reservoirs and its application for determining pressure variations and estimating production rate are carrying on for the past many years and are generally classified as numerical or analytical methods. Transport of shale gas in the reservoirs is a complex multi-scale transport process, which is from hydraulic fractures, i.e., macropores to the natural fractures, i.e., micropores (Javadpour et al. 2007a, b). Lots of researches have been done on transport mechanism of shale gas from matrix pores to the fractures. In general, most of the authors believe that the flow of gas in the fractures will follow Darcy’s Law, but the flow behavior of gas in matrix pores is still controversial. Zuber et al. (2002), Schepers et al. (2009), Wang and Reed (2009), Song and Ehlig-Economides (2011) and Song (Song and Yang 2013) conducted several studies and proposed that the flow of gas from the matrix pores to the fractures in shale gas reservoirs follows Darcy’s law.

Rushing et al. (1989), Dahaghi (2010) and Dahaghi and Mohaghesh (2011) have proposed that the flow of gas from the matrix pores to the fracture network is by diffusion. Javadpour (2009) and Ozkan et al. (2010) state that the flow and diffusion take place at the same time when the gas migrates from matrix pores to fracture network. As the permeability of the reservoir varies with location, it is not possible to have a unique permeability for the entire reservoir. For representing a uniform permeability for shale gas reservoirs, several investigations were performed on apparent gas permeability for representing the gas flow in shale reservoirs. Several investigations on apparent gas permeability have been done for representing the flow of gas in the nanopores (Clarkson and Nobakht 2011, Clarkson et al. 2012a, Clarkson and Williams 2012b; Michel et al. 2011; Civan et al. 2011; Sakhaee-Pour and Bryant 2012; Javadpour et al. 2007a, b; Javadpour 2009; Swami et al. 2012, Swami et al. 2013; Fathi et al. 2012).

One of the major factors in determining the productivity index of the shale gas reservoir depends upon the fracture network (Brown et al. 2009). In general, all the fractures are sourced by the matrix system. In most of the cases, a question arises about the contribution of shale matrix system to the fracture system. Unfortunately, with the available research a complete understanding of fluid transfer from shale matrix to fracture network is unknown. The present studies revealed that the main contributor to the flow of gas in the matrix is Darcy’s flow, which is induced due to pressure differential between the matrix and the fracture. Many authors have made different assumptions regarding the flow of gas in the shale matrix, as the fundamental assumption of Darcy’s flow in shale matrix revealed that the gas flow in the nanopores is considered negligible (Ozkan et al. 2010). In oder to have a clear idea about the flow of gas in the shale matrix, a detailed research has to be done.

Recently, Javadpour et al. 2007a, b; Javadpour 2009 described the flow in shale matrix by Knudsen diffusion and slip flow in nanopores, Darcy’s flow in the micropores, desorption from surface of the kerogen and the diffusion from the surface of the solid kerogen. Our objective in this paper is to include more detailed description of flow in shale matrix to the modeling of production from the fractured shale gas reservoir. Here, we limit our focus on Darcy’s flow, non-Darcy’s flow and desorption flow process. Desorption of gas in shale reservoirs has been linked to the coal-bed methane reservoirs where gas desorbs from the surface of the coal matrix block to the cleats (Induced Fractures). In shale gas reservoirs, the gas will be stored in the form of free gas and the adsorbed gas.

Here, we are presenting an updated dual-mechanism dual-porosity that accounts the free gas in the reservoir pores and the adsorbed gas on the surface of the kerogen. We consider a cubical matrix blocks, which consists of free gas and the adsorbed gas. As the pore space in the matrix reduces due to pressure reduction in the reservoir, the compressibility of the reservoir is also considered. The general formulation presented here represents the flow of gas in the matrix. The reservoir is divided into 5*5*5 matrix blocks. Now, mass balance equation is developed by considering a unique matrix block in the reservoir.


Shale gas is connected with significantly less carbon emissions as compared to coal. It can also decrease energy costs because huge amount of shale gas production would likely cause a decline in the price of natural gas. High shale gas production would also help our energy security and reduce our dependence on foreign fossil fuels (Ding et al. 2011). Shale gas could also provide better and cleaner energy option for many developing countries that are currently dependent on coal which is the dirtiest energy source.


There are additionally some disservices of shale gas. Shale gas, in spite of being essentially cleaner vitality source when contrasted with coal, regardless frees noteworthy carbon outflows, in this way being less satisfactory from ecological perspective than renewable wellsprings of vitality (Hong et al. 2013). Additionally, ecological danger as potential spillages of methane gas from different wells of shale gas could balance the decrease of carbon dioxide and atmosphere advantage of changing from coal to shale gas.

The fast improvement in shale gas businesses could back off the advancement of renewable vitality, particularly if shale gas gets to be one of the least expensive vitality choices accessible. Renewable vitality is thinking that it is hard to contend with coal, and with modest and effectively accessible shale gas, things could turn out to be much more terrible for the area of renewable vitality. Right now, the removing expense of shale gas is higher when contrasted with the expenses of extraction of routine gas or coal; however, the up-and-coming upgrades in boring innovations could diminish the extraction costs (Alahmadi 2010).


In the process of fluid flow characterization in shale reservoir, two basic approaches were used. The basic and initial approach is developing nonlinear partial differential equations which represent the flow of gas in the matrix and the flow of gas in the induced fractures and compiling these equations in MATLAB.

A second approach of solving and obtaining all the parameters will be used by the help of simulators. For the matter of credibility, the result of the equations derived from the first approach which are solved in MATLAB will be cross-checked with the results of the second approach using CMG-IMEX reservoir simulator.

Approach by MATLAB

MATLAB is used in our project to solve number of partial differential equations. The set of partial differential equations are solved by finite difference method by assuming some of the constants using the standard literature (Zhang and Yuan 2002). A generic equation is simplified which will change according to reservoir matrix in three dimensions by variables which are (i, j, k) which vary according to (x, y, z). The number of equations formed will depend on dimensions of the number of matrix assumed; for example, for n = 5, number of equations formed will be 5*5*5 = 125 equations. These equations are solved by using a MATLAB code using functions of matrices.

The following is the list of variables that are used in the code which can be later changed of different conditions:

  • Pm–For initial reservoir pressure.
  • T–For total number of days.
  • Dt–For time period.
  • dx, dy, dz–For reservoir length, breadth and depth.
  • N–For number of Matrix we want to solve.

A number of functions are created to facilitate the calculation of constants with respect to pressure changes at each reservoir point and with respect to time.

A nested loop is used to run the solution code by assigning the constants of each equation in a 3-D matrix and solving it for the values of the variables (Daniel Arthur and Coughlin 2012). A level 5 nesting codes are used in our coding. The final solution matrix is displayed using four-dimensional matrix for every time step.

Code description

The complete code is attached with Appendix 1.

The motive of this code is to solve a generalized linear equation for pressure values at each point in the given matrix. The number of unknown variables in the given matrix depends on the order of matrix assumed; for example, if we assume a matrix of the order of [5 × 5 × 5], then the number of elements in the given matrix will be 125, which further means that the number of unknown pressure points to be calculated by the generalized equation assumed previously would be 125.

The matrix that is considered in the project is in accordance with the dimensions of the physical shale rock matrix with following dimensions:

  • Length = 22ft.
  • Height = 12ft.
  • Thickness = 2ft.

To perform a solution for linear equations with such a large number of values, a generalized code is prepared which can be easily modified and used for different values for order of matrix, initial pressure and other different dependent variables.

The generic equation that was derived was a linear equation with seven unknown variables. A particular set of these unknown variables is unique for every point in the matrix. In this way by employing an equation for every point and calculating the corresponding seven unknown variables, pressure difference value at every point in the given matrix can be found for a particular value of time (Dreier 2004).

This step is repeated for every time interval for the complete duration of the project life or the well production period of the shale reservoir. The value of time interval and the complete duration of the project life are assumed as follows:

  • Time interval (dT): 10 days.
  • Time duration (T): 1000 days.
  • Number iterations done: 100.

The equation along with the seven unknown variables has their corresponding coefficients and a single constant value at the right-hand side of every equation. The value of these coefficients and constants depends on the pressure values of the matrix of the preceding time interval. These coefficients and constants change at each point in the matrix and with each time step. At the initial time, i.e., at t = 0, governing factor for these coefficients and constants is the initial pressure value. As the time changes, i.e., at the second time step the value of these will be depending on the previous pressure value at the respective point. So these coefficients and the constants are made dynamic whose values are getting updated with each successive iteration. A code snippet is attached (Fig. 1).

Fig. 1 Snippet showing the calculations of finite difference constants

Fig. 1 Snippet showing the calculations of finite difference constants.

Dynamic updating of coefficients and constants

The complete set of the equations, i.e., 125 equations, are solved using the standard matrix analogy. Every coefficient and constant of matrix are identified with the help of index position which corresponds to each point in the matrix as discussed above. All these equations are first arranged in a standard form and the left-hand side, i.e., the coefficient values are stored in a newly defined two-dimensional matrix, with every row containing the coefficients of that particular equation columnwise. The right-hand sides of the equations, i.e., the constants, are stored in the form of a single column matrix. This is done by using the following code (Fig. 2).

Fig. 2 Snippet showing the arrangement of equations in 2-D and column matrix.

Fig. 2 Snippet showing the arrangement of equations in 2-D and column matrix.

These set of matrices are then solved by the standard matrix form, i.e., AX = B. The inverse of the two-dimensional matrix is calculated and multiplied by the single column matrix to achieve 125 pressure values which are also in the form of a single column matrix. These single column matrixes with 125 fresh calculated values are then assigned to their respective places in the matrix using the technique of index assignment casting. This complete technique is shown in the following snippet (Fig. 3).

Fig. 3 Snippet showing the reverse allocation of pressure values to the original matrix.

Fig. 3 Snippet showing the reverse allocation of pressure values to the original matrix.

In the end, a 4-D matrix is considered with the fourth order to be made equal to the number of time steps, at which each set of 3-D matrices containing the pressure values is stored.

NOTE: The complete code of the project is attached with Appendices 1 and 2.

Approach by CMG-IMEX simulator

In this, we are going to present the approach for preparing the flow model by the help of validation software for simulation of CMG-IMEX simulator (Li 2007). This is a unique in its kind of software for showing the shale gas simulation at various points in the grid blocks.

The methodology basically involves of following steps to perform the analysis:

  1. I/O Control
  2. Reservoir
  3. Components
  4. Rock fluid
  5. Initial conditions
  6. Numerical
  7. Wells and recurrent


It is the second type of property set rather the main type of data set which is mainly composed of ten further parameters

  1. Grid: involves the following steps to perform the analysis grid type Cartesian–60*60*5 with a dual-porosity model, and the pinch out thickness of 0.0002 is set.
  2. Array properties (Figs. 4, 5, 6).   Open image in new window
  3. Rock fluid properties: default defined values.
  4. Sectors: default defined values.
  5. Aquifers: no aquifer is potentially used in various models.
  6. Lease plane: default defined values.
  7. Rock compressibility (Fig. 7).
  8. Compaction: default defined Values.
  9. Depletion: default defined Values.
  10. Flux sectors: default defined values.

Fig. 4   Array properties 1

Fig. 4   Array properties 1.

Fig. 5   Array properties 2

Fig. 5   Array properties 2.

Fig. 6   Array properties 3

Fig. 6   Array properties 3.

Fig. 7    Values input of rock compressibility

Fig. 7  Values input of rock compressibility

Emanuel Martin
Emanuel Martin is a Petroleum Engineer graduate from the Faculty of Engineering and a musician educate in the Arts Faculty at National University of Cuyo. In an independent way he’s researching about shale gas & tight oil and building this website to spread the scientist knowledge of the shale industry.

Leave a Reply

eleven − two =