Modeling and Simulation of Natural Gas Production from Unconventional Shale Reservoirs

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.


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.

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.
where u g 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 ϕ de- notes 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.

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.
where G is the potential releasable gas amount, P represents the average reservoir pressure, V L stands for Langmuir volume constant and P L 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.

(
) ( ) where p 0 represents the original gas pressure, v m denotes maximum adsorb gas volume for a single layer, and C symbolizes a heat constant which is defined in Equation ( 4).
where E 1 and E L 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.

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).
where 0.67 31.57 µ is the gas viscosity, c is the compressibility of gas, P reservoir pressure at current time step, k represents the permeability, and M g 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).
where k app 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.

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: 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).
In order to represent slip page flow in matrix blocks modified pseudo-pressure is given in Equation ( 11).
( ) For natural fracture pseudo-pressure: 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).

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 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.

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.

Summary and Conclusions
We presented our new semi-analytical reservoir model with incorporating fundamental flow mechanisms ob-   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:  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-linearity 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.

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

Figure 2 .
Figure 2. BET and Langmuir isotherm parameter analysis results.

Figure 3 .
Figure 3. Cumulative gas production with and without slippage flow.

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

Table 1 .
Reservoir parameters for history matching.

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

Table 3 .
Data used for parameter analysis.