WO 2013/089784 PCT/US2011/065566 SYSTEM AND METHOD FOR SIMULATION OF GAS DESORPTION IN A RESERVOIR USING A MULTI-POROSITY APPROACH Cross-Reference to Related Application(s) This application is a continuation-in-part, and claims priority to, PCT Application 5 No. PCT/US2011/44178, filed July 15, 2011, assigned to the assignee of the present application, and hereby incorporated by reference in its entirety. Background Reservoir simulation is an area of reservoir engineering that employs computer models to predict the transport of fluids, such as petroleum, water, and gas, within a 10 reservoir. Reservoir simulators are used by petroleum producers in determining how best to develop new fields, as well as in generating production forecasts on which investment decisions are based in connection with developed fields. Fractured reservoirs present special challenges for simulation because of the multiple porosity systems or structures that may be present in these types of reservoirs. 15 Fractured reservoirs are traditionally modeled by representing the porous media using two co-exiting pore systems or structures interconnected by flow networks, in what is referred to as dual porosity analysis. One type of pore system used in the prior art is the rock matrix, defined with matrix nodes, is characterized by high pore volume and low conductivity. The other type of pore system used in the prior art are induced fractures, and 20 defined with fracture nodes, is characterized by low pore volume and high conductivity. These prior art reservoir simulation methods and systems typically treat absorbed gas within the reservoir as residing in the rock matrix pores of the reservoir. For example, in one simulated representation, referred to as dual-porosity, single-permeability ("DPSP"), matrix simulation nodes communicate only with fracture simulation nodes, and the analysis 25 focuses on mass transfer and fluid flow of hydrocarbons between matrix nodes and fracture nodes. In DPSP, fracture nodes can also communicate with other fracures , which communicate with both matrix simulation nodes as well as other fracture simulation nodes. In another simulated representation, referred to as dual-porosity, dual-permeability ("DPDP"), matrix simulation nodes communicate with both fracture simulation nodes and 30 as well as other matrix simulation nodes, and the analysis focuses on mass transfer and fluid flow of hydrocarbons between matrix nodes and fracture nodes as well as between matrix nodes and other matrix nodes. 1 Those of ordinary skilled in the art will appreciate that "nodes" as used herein refer to the an elemental representation of pore structures within a simulated reservoir., while "zones" refer to a collection nodes within the simulated reservoir. Unknowns such as pressures and composition are solved for, typically on a node by 5 node basis, at desired time and/or depth increments. One particular type of reservoir encountered in oil and gas reservoir simulation is a shale reservoir. Shale reservoirs typically include large pores or vugs. Vugs are pore spaces that are comparatively larger than pore spaces of the rock 10 matrix. Kerogen resides in this system of vugs within the porous rock matrix. Vugs may or may not be connected to one another. "Separate vugs" are vugs that are interconnected only through the interparticle porosity, i.e., the rock matrix porosity, and are not interconnected to one another (as are matrix pore volumes and 15 fracture pore volumes). "Touching vugs" are vugs that are interconnected to one another. Because of their separate physical and mechanical characteristics, the fluid retention and transport properties of vug pore systems are different from those of both the matrix and fracture systems, and have not heretofore been adequately addressed with analysis utilizing only matrix porosity systems and induced porosity 20 systems. In other words, because of the geologic complexities of shale reservoirs, traditional dual porosity reservoir modeling techniques do not adequately predict mass transfer and fluid flow characteristics of shale reservoirs. SUMMARY OF THE INVENTION 25 In one aspect, the present invention provides a method for performing simulation of a shale reservoir represented by a model, the method comprising computer-implemented steps of: characterizing at least three different porosity types for the modeled reservoir, the at least three different porosity types selected from the 30 group consisting of natural fracture pore systems, matrix pore systems, induced fracture pore systems and vug pore systems; identifying transfer terms between the at 2 24/07/15,ag21678 amended speci pages,2 least three porosity types; and solving a linear system for the model using a linear solver. In another aspect, the present invention provides a computer program 5 product comprising non-transitory computer-readable medium having stored thereon instructions executable by a computer for causing the computer perform simulation of a reservoir represented by a model having a plurality of porosity nodes, the instructions for causing the computer to: characterizing at least three different porosity types for the modeled reservoir, the at least three different porosity types 10 selected from the group consisting of natural fracture pore systems, matrix pore systems, induced fracture pore systems and vug pore systems; identifying transfer terms between the nodes of at least three porosity types; and solving a linear system for the model using a linear solver. 15 In yet another aspect, the present invention provides a method for drilling one or more wellbores in shale reservoir, which method comprises: modeling an oil and gas shale reservoir having natural fracture pore systems, matrix pore systems, induced fracture pore systems and vug pore systems; characterizing at least three different porosity types in the modeled reservoir, wherein the characterized pore 20 types are utilized to create one or more subgrids that represent a zone within the reservoir, and wherein the zone includes at least one node of each porosity type; assigning transfer terms between the at least three porosity types, wherein transfer terms between nodes are assigned on a nodal basis to the nodes of the subgrids; and drilling a wellbore in accordance with the selected characteristic. 25 BRIEF DESCRIPTION OF THE DRAWINGS A more complete understanding of the present disclosure and advantages thereof may be acquired by referring to the following description taken in 30 conjunction with the accompanying figures, wherein: 2a 2 4 /0 7 /15,ag21678 amended speci pages,2 Fig. 1 illustrates an example of a reservoir simulation model comprising multiple wells. Fig. 2 illustrates a representation of an example formation comprising a 5 complex network of artificially-induced fractures. 2b 24/07/15,ag21678 amended speci pages,2 WO 2013/089784 PCT/US2011/065566 Fig. 3 illustrates a simulation grid of a formation comprising a highly deviated wellbore surrounded by natural fractures and a complex network of artificially-induced fractures. Detailed Description 5 To overcome the above-noted and other limitations of the current approaches, one or more embodiments described herein comprise a reservoir simulator including a unique manner of handling gas desorption in shale gas reservoir simulations by rigorously simulating the flow mechanism that occurs therein. It has been found that the mechanism for desorption of gas in a shale gas reservoir 10 is based on the existence of four separate porosity systems, each of which is incorporated in the method and system of the invention. In the method and system of the invention, each of these four porosity systems is separately characterized and incorporated into the model. The four porosity systems are the matrix porosity system, the induced fracture porosity systems, the natural fracture porosity system and the vug porosity system. As 15 explained above, heretofore, only the matrix porosity system and the induced fracture porosity systems have been used in reservoir modeling in the past. The method and system of the invention incorporate natural fracture porosity systems and vug porosity systems. The innermost of these porosity systems is the kerogen vugs, which contain the gas saturation as wetting fluid. The other porosity systems, which are the rock matrix, the 20 induced fracture network and the natural fracture network, function as conduits for the gas contained in the kerogen of the shale. Rather than residing in pores throughout the porous rock matrix, the adsorbed gas is generally found only in the kerogen vugs. Natural fractures exist near the vugs, which natural fractures may or may not be open. The framework rock matrix of the porous medium connects the complex natural fractures to the 25 hydraulic induced fractures near the well. The only accurate way to treat this system in accordance with the flow characteristics thereof, and in particular to account for the vugs throughout the matrix, requires a multiple porosity simulation system in which the vuggy portions of the formation containing the kerogen are connected to the rock matrix and the natural fracture system. The matrix and natural fractures are connected to the induced 30 fractures near the wellbore, which fractures will in turn have different properties than the natural fractures due to the presence of fracking fluid and perhaps proppant, and hence one 3 WO 2013/089784 PCT/US2011/065566 reason why natural fractures pore systems and induced fracture pore systems are characterized and separately analyzed in the method and system of the invention. The simulation system described in the aforementioned PCT Application No. PCT/US2011/44178 provides a unique tool for simulating general multi-porosity systems 5 in which fluid flow through several porosity systems is modeled using various equations and connectivities in accordance with characteristics of the porosity systems. The embodiments described herein utilize this feature to simulate in a unique fashion the desorption of gases from shale gas reservoirs. The equations provided hereinbelow are transfer functions derived from field observations and laboratory measurements of the 10 desorption process from the kerogen vugs, matrix and fractures of a reservoir. The transfer functions are then utilized by the simulation system to simulate the complex fracture system for the shale as reservoir coupled with the kerogen desorption. In one embodiment of the invention, vuggy porosity is used to model the kerogen desorption from within a complex fracture system comprised of both induced and natural fractures. 15 Fig. 1 is a block diagram of an exemplary computer system 100 adapted for implementing a reservoir simulation system as described herein. In one embodiment, the computer system 100 includes at least one processor 102, storage 104, I/O devices 106, and a display 108 interconnected via a system bus 109. Software instructions executable by the processor 102 for implementing a reservoir simulation system 110 in accordance with the 20 embodiments described herein, may be stored in storage 104. Although not explicitly shown in Fig. 1, it will be recognized that the computer system 100 may be connected to one or more public and/or private networks via appropriate network connections. It will also be recognized that the software instructions comprising the reservoir simulation system 110 may be loaded into storage 104 from a CD-ROM or other appropriate storage 25 media. In one embodiment, a portion of the reservoir simulation system 110 is implemented using reservoir simulation software. In this embodiment, a "subgrid" data type is used to offer a generalized formulation design. In one embodiment, this data type may be Fortran. The subgrid defines the grid domain and interconnectivity properties of 30 the nodes of the various porosity structures. It also tracks various node variables, such as pressure, composition, fluid saturation, etc. Subgrids are designated as being of a particular porosity type, e.g., natural fracture, matrix, induced fracture and vug. The nodes 4 WO 2013/089784 PCT/US2011/065566 that constitute these grids are correspondingly referred to as natural fracture nodes, matrix nodes, induced fracture nodes and vug nodes. Subgrids of different porosity types occupying the same physical space are said to be "associated". Connections between porosity types, and in particular, the nodes of the porosity types, are represented as external 5 connections, subgrid to associated subgrids. Internal (or intragrid) connections, and in particular, the nodes of a subgrid, represent flow connections within a porosity type. The modeling of a shale gas reservoir generally involves defining one or more elongated, highly deviated production wellbore, typically thousands of feet in length, with multiple hydraulic fracture zones disposed substantially perpendicular to the wellbore, 10 depending on the stress field in the formation. For certain formations, the stress field is such that a complex fracture system is induced between the large fractures emanating from the well. One representation of such fractures for an example formation is presented in Fig. 2 and designated by a reference numeral 200. The representation 200 has been derived from a finite element model of the porous media of the formation following high 15 pressure injection of fracture fluid and proppant. Heavier lines, such as those designated by reference numeral 202, represent fractures induced by hydraulic fracturing, as described above, and which have been modeled in the prior art, i.e., induced fracture porosity systems. Narrower lines and triangular features, such as those designated by reference numerals 204 and 206, respectively, represent a possible finite volume grid with which to 20 model the flow of fluid (primarily gas and water) in the complex fracture network and eventually to a horizontal production wellbore via the induced fracture, as illustrated in Fig. 3. Specifically, Fig. 3 illustrates a simulation grid 300 for an elongated, substantially horizontal wellbore 302, surrounded by induced fractures 304 and a complex natural fracture network 306. 25 It has been found that there are two features of the physics in reservoir simulation that are required to properly model shale gas flow to a wellbore in a reservoir: non-Darcy flow and gas desorption. In the method and system of the invention, these two features are considered when modeling mass transfer and fluid flow between the matrix nodes, natural fracture nodes, induced fracture nodes and vug nodes. 30 Non-Darcy flow is fluid flow that deviates from Darcy's assumption that fluid flow in the formation will be laminar. Non-Darcy flow is typically observed in high-velocity gas flow induced pressure differentials between the formation and the wellbore. 5 WO 2013/089784 PCT/US2011/065566 Specifically, when the flow at the wellbore reaches velocities in excess of the Reynolds number for Darcy (or laminar) flow, turbulent flow results and non-Darcy analysis must be utilized. The effect of non-Darcy flow is a rate-dependent skin effect. That is, as the velocity within the wellbore increases, there is an increase in the pressure drop between the 5 wellbore and the fracture. Thus, the typical equation for flow in the reservoir is modified to account for the effects of non-Darcy flow using the Forchheimer parameter P as shown in equation (1) below: (1) =( )q+Pp( )2 Ox KkA A 10 Where: P=pressure = pressure drop in a direction x ax p = viscosity K = permeability 15 kr = relative permeability A = cross sectional area to flow P = Forchheimer parameter p = density q = flow rate 20 For high velocity flow occurring in the fractures and near the wellbore, non-Darcy flow results in a significant increase in the pressure drop and therefore plays an important role in properly modeling shale gas production. Because prior art techniques did not tend to model natural fracture pore systems, to the extent non-Darcy flow analysis has been used in the past for reservoir modeling, it has only been utilized to model flow in matrix 25 pore systems and induced fracture systems. Unfortunately, inclusion of the effect of equation (1) in reservoir fluid flow represents significantly more effort than was required for the skin factor. Since velocity depends not only on pressure drop but also on viscosity and relative permeability, a highly non-linear dependence is added to the flow equations for gridblock-to-gridblock or 30 fracture-to-fracture non-Darcy flow treatment. The skin factor only requires a minor modification to the coefficient for the pressure loss between the wellbore and the reservoir 6 WO 2013/089784 PCT/US2011/065566 or fractures. Inclusion of the non-Darcy effect adds a significant non-linear term to the pressure equations and requires that this term be included in the linearization for the Newton-Raphson iteration to solve for the flow in the wellbore and reservoir. In turn, this may increase the number of non-linear iterations and therefore increase overall 5 computation time for the reservoir simulation. Gas desorption for shale development is an important, heretofore underutilized parameter in shale formation modeling. It is estimated that in some shale formations, more than 50% of the gas production will be due to desorption. To the extent gas desorption has been modeled in shale reservoirs, its use has been limited to desorption from the shale 10 matrix. It has not heretofore been applied to desorption analysis from kerogen vugs. Because economics are highly dependent on ultimate recovery from the formation, gas desorption must be treated for a shale gas reservoir simulation system to have any credibility. Moreover, it must be applied in a way that accounts for the existence of kerogen vugs within the reservoir. Desorption is described by the Langmuir equation 15 (equation (2) below) for isothermal desorption characteristics: V, .P (2) V =__ * PL+ P Where: Vg volume of gas contained in the porous medium VL = asymptotic adsorption volume 20 PL = pressure at which the adsorbed volume reaches VL P = reservoir pressure Use of equation (2) in the simulator results in a modification similar to that for dual porosity, single permeability ("DPSP"), in which a source of gas, i.e., vug nodes, are included in each grid, the volume of which depends on the change in matrix pressure over 25 a timestep in the simulator. For more rigorous treatment of the physics, consideration of sorption time and desorption effects on formation permeability may also be necessary. Sorption time is the time it takes for 63.2% of the gas to be desorbed as calculated using equation (2). In the case of shale gas, this time is generally extremely short and can be ignored. Similarly, the 30 effect of desorption on matrix permeability is generally very small for shale gas and can also be easily ignored. 7 WO 2013/089784 PCT/US2011/065566 In practice, when pressure is lowered in the horizontal production wellbore 302 (Fig. 3), the pressure is almost instantaneously lowered in all of the fracture system (including both induced fracture and natural fractures) to which the wellbore is connected. For those fractures connected directly to a kerogen vug, the pressure is also reduced from 5 the initial pressure. From equation (2) above, the kerogen must release gas into the surrounding reservoir fractures and matrix based on the Langmuir equation. Although the VL and PL parameters can be experimentally determined, often these are estimated by analogy and rules of thumb. Flow through the multiple fractures and matrix represents a significant difference with conventional treatment in which only the matrix contains the 10 adsorbed gas which is directly in contact with fractures. The more complex treatment using vugs should allow a more realistic simulation of desorption than is currently achieved since the complex geometry of the porous medium is correctly characterized. With reference to Figure 4, a flowchart is shown illustrating the steps of the process of the invention. The process is utilized to model flow characteristics to a wellbore of a 15 shale reservoir having kerogen vugs and is preferably performed in conjunction with a three dimensional model of a reservoir. In step 400, reservoir characterization is undertake in which at least three, and preferably four different pore types are described based on fractured shale characteristics. In one embodiment, at least three different pore types are identified, selected from the group consisting of natural fracture pore systems, matrix pore 20 systems, induced fracture pore systems and vug pore systems. In one embodiment, four different pore types are identified, namely natural fracture pore systems, matrix pore systems, induced fracture pore systems and vug pore systems. In any event, in step 402, the pore types are utilized to create one or more subgrids that represent a zone within the reservoir. Each zone includes a plurality of nodes of at least one of the pore types. In one 25 embodiment, a subgrid for at least three different pore types is created for a zone. In one embodiment, a subgrid for each of the four pore types is created for a zone. In step 404, once the shale reservoir has been described sufficiently, connectivities or transfer terms, if any, between the nodes are identified and assigned. This may include connectivity between similar nodes within the same subgrid, such as between matrix nodes 30 within a subgrid, or may include connectivity between the nodes of one subgrid and the nodes of another associated subgrid, such as between vug nodes and natural fracture nodes or between matrix nodes and vug nodes. These transfer terms are the parameters that effect flow rates among the various porosity types, such as, for example, initial pore pressures, 8 WO 2013/089784 PCT/US2011/065566 basis to the nodes of the subgrids. In one embodiment, the model consists of at least three different pore types and associated volumes which contain fluids which are to be modeled. In one embodiment, the model consists of at least four different pore types and associated volumes which contain fluids which are to be modeled. 5 In step 406, know magnitudes for the transfer terms may be assigned, such as, for example, densities, volumes, flow rates and compressibilities. In step 408, source terms are now incorporated as boundary conditions to the model in such a way that extraction of the gas is consistent with the wellbore's induced fractures. Put another way, to initiate flow in the simulation, a wellbore pressure is selected and 10 incorporated into the model. This pressure affects the flow in the induced fractures, which, in turn by virtue of the transfer terms, affects flow between the other porosity types. In step 410, a linear solver is utilized to solve for any unknown magnitudes of the transfer terms associated with the nodes. In one embodiment, non-linear equations are selected to model the reservoir and the subgrids and nodes thereof. In one embodiment, 15 the linear solver methodology is applied by subgrid or associated subgrids. The Newton Raphson method is then applied to linearize these non-linear equations. The linear solver then can be applied to the linear equations to solve for the unknowns. In one embodiment, this step may be iterated utilizing the resultant magnitudes until a desired degree of convergence is achieved between the linear and non-linear equations. 20 In step 412, optionally, once a desired degree of convergence is obtained and the magnitudes of the unknowns are identified, time may be incremented and/or the wellbore parameters, such as the boundary conditions of pressure, may be altered to achieve a desired level of mass transfer and fluid flow for the modeled reservoir. The foregoing methods and systems described herein are particularly useful in 25 drilling wellbores in shale reservoirs. First a shale reservoir is modeled as described herein to design a well completion plan for a well. In an embodiment, the drilling well completion plan includes the selection of a fracturing plan, which may include the selection of fracture zones and their positioning, fracturing fluids, proppants and fracturing pressures. In other embodiments, the drilling well completion plan may include selecting a 30 particular trajeocry of the wellbore or selecting a desired wellbore pressure to facilitate mass transfer and fluid flow to the wellbore. Based on the model, a drilling plan may be 9 implemented and a wellbore drilled in accordance with the plan. Thereafter, in one embodiment, fracturing may be carried out in accordance with the model to enhance flow from the reservoir to the wellbore. In another embodiment, wellbore pressure may be adjusted in accordance with the model to achieve a desired degree of mass 5 transfer and fluid flow. Those of ordinary skilled in the art will appreciate that while the method of the invention has been described statically as part of implementation of a drilling plan, the method can also be implemented dynamically. Thus, a drilling plan may be implemented and data from the drilling process, and in particular, the actual flow characteristics of the reservoir, may be used to update the model for the 10 drilling of additional wellbores within the reservoir. After implementing the drilling plan, the system of the invention may be utilized during the drilling process on the fly or iteratively to calculate and re-calculate connectivity characteristics of the reservoir over a period of time as parameters change or are clarified or adjusted. In either case, the results of the dynamic calculations may be utilized to alter a 15 previously implemented drilling plan. For example, the dynamic calculations may result in the utilization of a heavier or lighter fracturing fluids. While certain features and embodiments of the invention have been described in detail herein, it will be readily understood that the invention 20 encompasses all modifications and enhancements within the scope and spirit of the following claims. Furthermore, no limitations are intended in the details of construction or design herein shown, other than as described in the claims below. Moreover, those skilled in the art will appreciate that description of various components as being oriented vertically or horizontally are not intended as 25 limitations, but are provided for the convenience of describing the invention. It is therefore evident that the particular illustrative embodiments disclosed above may be altered or modified and all such variations are considered within the scope and spirit of the present invention. Also, the terms in the claims have their 30 plain, ordinary meaning unless otherwise explicitly and clearly defined by the patentee. 10 24/0 7 /15,ag21678 amended speci pages,10 Throughout this specification and the claims which follow, unless the context requires otherwise, the word "comprise", and variations such as "comprises" and "comprising", will be understood to imply the inclusion of a stated integer or step or group of integers or steps but not the exclusion of any other integer or step or 5 group of integers or steps. The reference to any prior art in this specification is not, and should not be taken as, an acknowledgment or any form of suggestion that the prior art forms part of the common general knowledge. 10 10a 24/07/15,ag21678 amended speci pages,10