A review on the recent advances of the flow and transport phenomena in tight and shale formations is presented in this work. Exploration of oil and gas in resources that were once considered inaccessible opened the door to highlight interesting phenomena that require attention and understanding. The length scales associated with transport phenomena in tight and shale formations are rich. From nanoscale phenomena to field-scale applications, a unified frame that is able to encounter the varieties of phenomena associated with each scale may not be possible. Each scale has its own tools and limitations that may not, probably, be suitable at other scales. Multiscale algorithms that effectively couple simulations among various scales of porous media are therefore important. In this article, a review of the different length scales and the tools associated with each scale is introduced. Highlights on the different phenomena pertinent to each scale are summarized. Furthermore, the governing equations describing flow and transport phenomena at different scales are investigated. In addition, methods to solve these equations using numerical techniques are introduced. Cross-scale analysis and derivation of linear and nonlinear Darcy’s scale laws from pore-scale governing equations are described. Phenomena occurring at molecular scales and their thermodynamics are discussed. Flow slippage at the nanosize pores and its upscaling to Darcy’s scale are highlighted. Pore network models are discussed as a viable tool to estimate macroscopic parameters that are otherwise difficult to measure. Then, the environmental aspects associated with the different technologies used in stimulating the gas stored in tight and shale formations are briefly discussed.

1. Introduction

The problems associated with the scarcity of energy resources are intensified by the increased level of demands. Human beings are nowadays consuming more energy resources than ever before. Three fossil fuels, namely, petroleum, natural gas, and coal, have provided more than 80% of total US energy consumption for more than a century, EIA [1]. The depletion of energy resources from fossil fuels has reached alarming levels that necessitate decision-makers, research institutes, and industry to search for alternative energy resources. It seems, however, that mankind is not yet in the position to abandon totally his dependence on energy from fossil resources. In 2015, the renewable share of energy consumption in the United States was at its largest at nearly 10%, EIA [1]. The vast majority of our energy demands are satisfied by energy from fossil resources. Other energy resources (e.g., from renewable sources) may not be sufficient to supply our energy needs. Furthermore, most of our machinery, transportation, and devices are adapted to utilize energy from fossil resources. Therefore, there is still a trend to continue draining resources of fossil fuels to the last drop. Oil and gas resources can be divided in terms of their accessibility into conventional and unconventional resources. In the last decade, when the price of the oil exceeded considerably $100, it became economic to search for oil/gas resources in hardly accessible reservoirs using unconventional technologies. Although many of these projects have, nowadays, stopped due to the current decline in oil price, they are expected to resume once the price of the oil climbs up again. Based on projections in the US production of oil and gas from unconventional resources, tight oil production is expected to reach 7 million barrels per day and shale-gas production is expected to reach 79 billion cubic feet per day in 2040. Such increase in the production of tight oil and shale gas is driven by technological improvements that have reduced drilling costs and improved drilling efficiency. Figure 1 shows projected growth of the US shale-gas production until 2040, EIA [2]. New technologies have indeed increased our capacity to acquire resources that have once been thought inaccessible. In other words, energy resources that are easily accessible no longer exist and, nowadays, the search has been for resources that require larger investment and new technologies that are able to mine and to transport these resources at reasonable costs.

In the oil and gas industry, the types of oil and gas deposits are classified into what are called plays (Figure 2). Such plays are categorized based on many factors including the geology and the technology required to produce the oil. As seen in Figure 2, the fringe regions surrounding the area of historical production (halo zone) are likely to contain oil. The reason they have not been included with the early production region may be that the geologic properties are not as favorable as those within previously producing areas. Different technologies are used for different plays. To displace the oil contained in the halo zones, new technologies including horizontal wells have been used. The geologic formations in the halo zone are generally characterized by lower permeability. Therefore, oil reservoirs, in which geologic formations are of lower permeability (i.e., tight), require newer technologies to get access to oil reserves and also to be able to displace the oil. The distinction between conventional and unconventional reservoirs has been, in most cases, based on formation permeability. Unconventional reservoirs are acquired using horizontal wells as compared with vertical wells that have been traditionally used in conventional reservoirs. The purpose of drilling a horizontal well as compared with vertical wells is to increase the contact area between the reservoir and the wellbore. The horizontal leg of the well can extend up to 5 km. To stimulate tight oil reservoirs once the well has been drilled, hydraulic fracking is used. In this process, fluids are pumped into the wellbore at very high pressure to open existing fractures or to create new ones. Through these fractures, oil can flow to the wellbore (Figure 3). The types of fracture fluids vary depending on the reservoir formations with the water representing the base fluid. Additives (0.5% to 2% of the total fracturing fluid volume) are added to the water to reduce friction, control microorganism growth, and prevent corrosion. To maintain the fractures open during the production, sand particles (proppants) are pumped with the fracking fluids. The volume of fracking fluids and the amount of proppant used of hydraulic fracking vary depending on the required rate of production.

2. Characteristics of Tight and Shale Formations

Tight and shale formations are rocks with pores so small or poorly connected that the oil and natural gas cannot flow through them easily. In these formations, hydrocarbons in the form of crude oil, natural gas, and natural gas liquids may exist in considerable quantities. Like all hydrocarbons, they formed over millions of years, when organic material (e.g., plants and microorganisms) was buried and subjected to increasing heat and pressure and slowly transformed to oil and natural gas. Some of these hydrocarbons escaped into adjacent rock layers that are relatively easy to extract because of their relatively higher porosity and permeability. However, the majority remained locked in tighter, lower permeability layers where they could not be extracted through conventional means. The classification of petroleum reservoirs into conventional and unconventional is in part related to the technology used to extract entrapped oil and gas. Conventional oil reservoirs are those where the geologic formations are characterized by relatively higher permeability that allows easy transport of oil/gas reserves. On the contrary, unconventional reservoirs are characterized by permeability which is much lower such that conventional techniques to displace the oil/gas may not work. Figure 4 shows a schematic diagram of the classification of geologic formations according to their permeability.

Shale formations are one of the tightest rock formations in terms of their permeability. Most of the oil and gas reserves exist in an organic portion of rock mass called kerogen. Figure 5 shows an aerial view of shale sample where a kerogen area exists. The kerogen region is highly porous where shale gas exists in the pore space and adsorbed at the surface. Organic matter consists of kerogen (~90%) and bitumen (~10%). Figure 6 shows a schematic diagram for the composition of typical shale formation. As indicated by Bohacs et al. [3], Prasad et al. [4], Passey et al. [5], Wang and Cao [6], and others, when kerogen matures, it produces oil and gas. Kerogen does not have a specific structure as it is a mixture of organic materials in which the chemical compounds can vary from one sample to another. Kerogen maturation occurs when it is subjected to higher temperature for longer periods of time. Thermal decomposition breaks small molecules leaving behind a more resistant kerogen residue. The smaller molecules become eventually natural gas.

3. Characteristic Length Scales

Porous media exist almost everywhere around us from naturally occurring media to manmade systems. They encounter a larger spectrum of length scales ranging between global and regional scales all the way towards micro- and even nanoscales. Groundwater, petroleum, and geothermal reservoirs are examples of such larger size domains whereas cells, membranes, and living organisms are examples of micro- and nanosize domains. Porous media applications even extend beyond our planet towards nearby terrestrial planets to explore, for example, the existence of water underneath the surface in geologic formations. Porous media applications also span quite a large spectrum of time scales, from phenomena that take very long time scales to be of noticeable influence (e.g., groundwater flows) to others that take quite shorter period of time (e.g., flows in fluidized bed reactors). Therefore, no wonder there is extensive interest devoted to characterizing, understanding, and modeling several phenomena in porous media. Having such a wide spectrum of length and time scales, there exist several frameworks to handle such rich systems. One may be able to highlight three such frameworks. These are molecular scale simulations, pore-scale simulations, and continuum scale simulations. In all these simulation methodologies, upscaling techniques are required to produce integral variables that can easily be determined and measured. Two length scales are important in determining which framework may be used. These are the length scale characteristic to the domain and the length scale characteristic to the heterogeneity. Therefore, if the length scale characterizing the domain, L, is much larger than that associated with the scale of heterogeneity, d, (i.e., if ), the continuum hypothesis may apply (Whitaker [79]; Gray [10]; Hassanizadeh and Gray [11]; Cushman [12]; Carbonell and Whitaker [13]; Bachmat and Bear [14]; Quintard and Whitaker [15]). If, on the other hand, the length scale associated with the domain is on the order of the scale of heterogeneity (i.e., if ), the continuum approach may not be appropriate and a more detailed description at the pore scale may be required. In some situations, it so happens that the characteristic length scale is on the order of the average distance between molecules of the fluids such that even pore-scale simulation may not be correct. In this case, molecular simulations may be the appropriate framework. Figure 7 shows a schematic representation of the various length scales associated with porous media. The first graph to the left is an example of a domain that can be treated as a continuum, the middle graph may be treated using pore-scale modeling approaches including pore network models, and the last graph may be studied using molecular simulation techniques. To establish the continuum hypothesis, an averaging volume (called representative elementary volume, REV) over which upscaled quantities are determined needs to be defined. Such averaging volume is chosen such that upscaled quantities are free from scaling variations. To establish such requirements, the length scale characterizing the REV, , should be large enough compared with pore-scale heterogeneity, d, and small enough compared with the size of the domain, L (Salama and Van Geel [16, 17]). Therefore, Several approaches have been used to derive the equations governing flow and transport in porous media. These include the method of volume averaging, theory of homogenization, and theory of mixtures. Within the method of volume averaging, the pioneering works of Whitaker and his group, Bear and his group, Gray and Hassanizadeh, and others have paved the road for the advancement in the study of flow and transport in porous media. Salama and Van Geel highlighted the notion that the averaging process can be understood by requiring that the amount of any conservative quantity within any control volume (larger than or equal to the REV) must be the same whether calculated using pore-scale equations or upscaled one. To facilitate the analysis, Hassanizadeh and Gray [11] introduced the phase function which is defined such that in the -phase and elsewhere. Salama and Van Geel [16, 17] postulated that if is an arbitrary volume such that , then the following relationships apply:where is the density of the -phase, is any intensive property per unit mass (e.g., energy and momentum per unit mass), is the volume fraction of the -phase, and the quantities in are averaged over the volume of the -phase. When the volume is the REV, such integrations produce formulas for upscaling. When the characteristic length of the domain of interest is on the order of the average distances between molecules (sometimes the mean free path is used instead), then the collision of molecules with the boundaries becomes significant. In other words, in such cases, the molecules spend more time in the vicinity of the walls rather than in the bulk. Under these circumstances, the assumption of thermodynamic equilibrium becomes questionable. That is, momentum and energy transport and the convergence to equilibrium are based on the collisions between molecules in the bulk fluid, which no longer exist at higher Knudsen number. In such small size domains, even the definitions of macroscopic variables (e.g., density, pressure, and temperature) as a manifestation of the average of the behavior of fluid particles within a representative volume may not be unique and will essentially be size-dependent. Under these conditions, the tools of molecular simulation become more appropriate to handle the state of equilibrium rather than classical bulk-phase thermodynamics. The Knudsen number is used to characterize when the continuum approach fails to apply in fluids. Knudsen number is the ratio of the mean free path and the characteristic length scale of the study domain (i.e., ). Figure 8 shows a map of the different flow categories when the length scale of the domain decreases (i.e., increasing Knudsen number).

4. Thermodynamics Associated with Transport Phenomena in Tight Formations

Classical thermodynamics is based on a set of postulates that determines the state of the system under different conditions. If a system is disturbed, its constituents interact in an attempt to return to its initial state or to establish a new equilibrium state. The constituents of the system interact with each other and across the boundary with the surroundings in a definite manner that is determined by the laws of thermodynamics. The collisions of the molecules of the system in the bulk help homogenize local disturbances quickly. When the characteristic length scale of the system is on the order of the average distance between the particles, the collision with the walls of the system becomes dominant. Under these circumstances, a number of interesting phenomena have been observed: () momentum transfer is increasingly controlled by the wall collisions, () the fluid and flow properties start to fluctuate in a selected differential volume due to the lack of a sufficient number of molecules needed for statistical accuracy, () gas attains finite velocities in the close proximity of the wall (i.e., gas slips over the wall), and () thermodynamic variables like pressure drop, shear stress, heat flux, and corresponding mass flow rate cannot be predicted from flow and heat transfer models based on the continuum hypothesis. The appropriate flow and heat transfer models depend on the range of the Knudsen number as depicted in Figure 8. These and others necessitate that all the assumptions of classical thermodynamics need to be revisited.

5. Pore-Scale Phenomena

Natural porous materials abound in the earth, such as expansive bentonite clays and porous limestone, and in living organisms, such as the mastoid bone with its porous air cavities. Extraction of natural gas from pores in shale by hydraulic fracturing has transformed the energy agenda of the US and the whole world (Kobek et al., 2015). Understanding the properties of these pores in terms of their structure and transport properties is now emerging as a grand challenge in terms of the scaling from the microscopic or nanoscale regime to the macroscopic world. This transition from the nanoscale to a macroscopic world is known as mesoscale science, in which the field of pore-scale phenomena is now emerging as one of the frontiers of science and many engineering disciplines and dominating the extraction of these resources. Permeability is one of the most fundamental properties of any reservoir rock required for modeling hydrocarbon production. However, shale permeability has not yet been understood fully because of the complexities involved in modeling flow through pore throats (Sakhaee-Pour et al., 2012). New pore-scale models with a reservoir simulation algorithm to predict gas production in gas-bearing shales have been proposed in recent years, especially concentrating on the permeability determination with the development of shale-gas engineering, which simultaneously consider the effects of slip flow (Klinkenberg effects) and Knudsen diffusion [18, 19].

6. Klinkenberg Effect in Shale

Slippage of fluid continuum upon encountering a solid surface has been investigated by Navier since the middle of the nineteenth century when he proposed an extended length where the velocity profile extrapolates to zero. Such extended length is, generally, so small in regular bulk flows in which the characteristic length scale is considerably larger than the average distance between particles. In micro- and nanochannels, however, such condition may not be satisfied and flow slippage may become pronounced. In shale formations, flow slippage in nanosize pores and fractures needs to be accounted for. Due to the gas-slippage effect, the permeability of a sample to a gas varies with the molecular weight of the gas and the applied pressure, which was first proposed by Klinkenberg [20] and so called Klinkenberg effect thereafter. He determined that the slippage of gases along the pore walls gives rise to an apparent dependence of permeability on pressure, which could be concluded as that liquid permeability () is related to gas permeability () bywhere is the mean flowing pressure and is a constant for a particular gas in a given rock type. This non-Darcy effect occurs when the mean free path length of the gas molecules is close to the average size of pores in a porous medium. This condition results in the acceleration of individual gas molecules along the flow path [21]. The Klinkenberg effect is especially important in low-permeable rocks, so it attracted more and more attention with the development of shale gas in recent years (Civan [22], Tanikawa and Shimamoto [23]). Wu et al. [24] proposed a set of new analytical solutions developed for analyzing steady-state and transient gas flow through porous media including Klinkenberg effects, which have been used to design new laboratory and field testing techniques to determine the Klinkenberg parameters. Other popular approaches for determination of shale-gas permeability are presented by Pazos et al. [25] and Sakhaee-Pour et al. (2012).

Experiments have been designed and performed to measure the shale permeability, to verify the effect theory, and to obtain the Klinkenberg parameter, which can be found in Jones [26], Faulkner and Rutter [27], Tanikawa and Shimamoto [28], Cui et al. [29], and Davarzani et al. [30]. Generally, the pore network is constructed based on scanning electron microscope (SEM) images and a drainage experiment in shale. Currently, there are three methods in common use for determining permeability of very low permeability rocks in the laboratory. These include (i) studying the permeation of inert gas (e.g., helium) through a core sample under either falling pressure or steady-state pressure techniques, (ii) building a digital realization of the core sample and performing a CFD analysis to determine the overall resistance, and (iii) using mercury (Hg) intrusion curves (from Hg porosimetry). Permeability measurements on rock core under confined conditions have been routinely used for conventional oil and gas reservoirs for over 50 years. However, traditional steady-flow permeability measurement (American Petroleum Institute (API) 1998) on core samples of very tight rocks such as most gas shales and coal is not practical because of the time scales involved and the instrumentation requirements for measuring extremely small pressure drops or flow rates. Another method for approximating permeability is by using Hg injection curves from Hg porosimetry. The relationship between Hg injection curves and permeability has been investigated by a number of authors (Thomeer, 1960, 1983; Swanson, 1981; Kamath, 1992, Carles et al., 2007). In Swanson’s (1981) method, for example, permeability is calculated by considering the Hg saturation and capillary pressure at the apex of a hyperbolic log–log Hg injection plot. Swanson (1981) developed and calibrated the relationship between permeability and Hg intrusion data from a suite of sandstone and carbonate samples. Hg intrusion as a permeability tool is not further considered in this paper, except to note that permeability or diffusion measurements on unconfined samples are at best instructive because permeability is known to vary with effective stress by several orders of magnitude (i.e., Bustin, 1997, [31]). Another traditional technique of measuring shale permeability using crushed samples was designed by Luffel and others (1992 and 1993) to measure matrix permeability only by eliminating natural and drilling induced microfractures. Although drilling induced fractures are common (Boyer and others, 2006), they can be minimized by selecting core plugs at locations without drilling induced fractures. Because pore networks in organic matter are most likely connected through microfractures, the connectivity of organic pore network can be significantly reduced in crushed samples. Although pore networks in organic matter and natural microfractures are important properties of shale and critical to shale-gas production, they are too small to be properly quantified in the laboratory or reservoir simulation. The compromissory but easy way is to include them in core permeability measurements as part of a lumped permeability value. By including important organic and microfracture pore networks, this lumped permeability can characterize gas shales better than the true matrix permeability. Theoretically enhanced models incorporating the Klinkenberg effect and using effective method to calculate the permeability of gas transport in gas-bearing shale formations have been presented to govern gas flow in shale reservoir, as shown by Civan [32, 33], Civan et al. [34], and Al-Bulushi et al. [35].

7. Knudsen Diffusion in Shale

When the mean free path of gas molecules is on the same order as the tube dimensions (as in shale), Knudsen diffusion, which is a typical kind of free-molecule diffusion, becomes important. Due to the influence of walls, Knudsen diffusion includes the effect of the porous medium. Gas mass flux by diffusion with negligible viscous effects in a nanopore is described as [36]where is molar mass, is the Knudsen diffusion constant, (=8.314 J/mol/K) is the gas constant, and is absolute temperature in Kelvin. The Knudsen diffusion constant is defined as [37]where represents the molecular weights of gas and is the mean pore size of the porous media. Many experiments have been performed to measure the Knudsen diffusion constant, which can be found in (exclusive for those listed in the above session) Reinecke and Sleep [38], Jarvie [39], and Freeman et al. [40]. Different algorithms to simulate gas diffusion have also been developed in [41, 42]. Malek and Coppens [43] studied the effects of surface roughness on Knudsen regime diffusion in porous media. Knudsen diffusion is a result of collisions of gas molecules with the pore walls, rather than intramolecular collisions, so we always consider that Knudsen diffusion and molecular diffusion compete with one another by a “resistances in series” approach [44]. Welty et al. [45] presented the different types of diffusion as shown in Figure 9. Combining Darcy’s law in unconventional systems with the total mass flux formula in the shale, we obtain the formula for the apparent permeability, which replaces the intrinsic permeability in conventional media. A commonly used formula for computing the apparent permeability of shale based on Klinkenberg effects and Knudsen diffusion is [46] where and . Here, is the intrinsic permeability (i.e., the permeability for sufficiently large pressure, sometimes called liquid permeability), is the apparent permeability, is the average pore size diameter, τ is the effective tortuosity of the pores (if the pore spaces are assumed to be straight and cylindrical capillaries, the effective tortuosity is equal to 1), is the universal gas constant, with the unit of J/mol/K, M is molecular weight of the fluid, with the unit of kg/kmol, α is the tangential momentum accommodation coefficient, dimensionless, is the dynamic viscosity of the fluid, and is the temperature in K.

We would like to comment that the above Knudsen diffusion description fails to capture systems having very small pores on the order of 1 nm such as in mature kerogen. In this situation, a gas molecule has a size similar to the pore dimension. This situation is very different from the Knudsen regime. In fact, in this situation, conventional fluid viscosity no longer makes sense because gas molecules interact essentially with the pore walls instead of interacting with gas molecules themselves. When the size of the gas molecules is only slightly smaller than the pore dimension, the diffusion rate that occurred is usually much larger than predicted by Knudsen diffusion because of the superlubricity effect. When the size of the gas molecules becomes larger than the pore dimension, the diffusion rate is rapidly reduced to zero because of the molecular sieving effect.

8. Pore Network Models

As indicated earlier, in tight and shale formations, the permeability is quite small such that reliable measurements may not be readily possible. Although the flow in a bulk sample of the shale rock is difficult to measure, it can be calculated if the internal structure of a rock sample is known. With the advancement in imaging technologies, nowadays, it became possible to construct a realization of the internal structure of rock samples using series of images. X-ray computed microtomography (micro-CT) is used to produce slice pictures of the rock sample that are then used to construct a pore network model mimicking the real pore structure (Figure 10). The flow in these systems can easily be studied using, for example, the Hagen-Poiseuille approximation (for single-phase flows), which can then be used to estimate the absolute permeability of the rock sample. Furthermore, the study of the problem of drainage and imbibition can also be considered using the pore network to estimate the relative permeability characteristics of the rock sample.

Several experimental techniques have been explored to provide three-dimensional details of the microstructures of rock samples. Computed tomography is a nondestructive imaging technique used to characterize the internal structure of several things including rock samples. Three types of CT systems are in common use, namely, medical CT, industrial X-ray generation tube, and synchrotron microtomography (Hazlett [47], Wildenschild et al. [48], Withers [49], and Schlüter et al. [50]). The typical spatial resolution range that medical CT systems can achieve may be between 200 and 500 microns, industrial X-ray tube systems range from 50 to 100 microns, and synchrotron based systems are from 1 to 50 microns. A recent review on these imaging techniques can be found in Blunt et al. [51]. On the other hand, other techniques including focused ion beams (FIB) and scanning electron microscopy (SEM) are essentially destructive (Tomutsa et al. [52], Curtis et al. [53], and Lemmens et al. [54]). SEM can be used to extract two-dimensional planar images of the microstructures. However, they do not provide the third spatial extent of the sample which is essential to determine connectivity regions. While FIB is very effective in generating higher resolution three-dimensional images, it is very much time-consuming due to the refocusing and repositioning requirements. A combination of FIB and SEM is usually used to compute structural details of porous media (Michael et al. [55], Keller et al. [56]). Nuclear magnetic resonance imaging (NMR), on the other hand, allows the imaging of the interior of the rocks to obtain the spatial distribution across a much larger scale (Callaghan [57], Blümich et al. [58]). NMR has the advantage that it requires shorter measurement time compared with other methods, which allows the analysis of larger quantities of samples to characterize field-scale hydraulic properties. For porosity measurements, mercury intrusion technique is probably the most popular for characterizing porous materials with pore sizes ranging from 3 nm to 500 μm (Giesche [59], Léon and León [60], and Rouquerol et al. [61]). Gas adsorption, on the other hand, is based on the adsorption behavior of the porous material, which, in turn, is a function of microstructural characteristics of the porous material. Traditionally, nitrogen (N2), argon (Ar), and carbon dioxide (CO2) are frequently used as adsorbates depending on the nature of the porous materials (Ravikovitch et al. [62], Groen et al. [63], and Settano et al. (2009)). A recent review on the methods discussed earlier can be found in Xiong et al. [64]. The next step after preparing the stack of images that describe the spatial distribution of the real rock sample is to build a three-dimensional realization of the pore space in the network. As reported by Xiong et al. [64], there exist a number of techniques to construct a pore network model (PNM) representing a given porous medium. They may be categorized into three methods, namely, statistical reconstruction (Adler and Thovert [65], Levitz [66], Roberts and Torquato [67], Ioannidis and Chatzis [68], and Manwart et al. [69]), grain-based model (Bryant et al. [70], Bakke and Oren [71], Lerdahl et al. [72], and Oren and Bakke [73]), and direct mapping model (Al-Raoush and Wilson [74], Jiang et al. [75], Shin et al. [76], and Raoof and Hassanizadeh [77]). Once the pore network structure has been established, different interesting phenomena can be investigated. These include the estimation of absolute permeability, relative permeability, adsorption, dissolution, and precipitation, biomass growth, and others. In shale and tight formations, the absolute permeability is estimated using the Hagen-Poiseuille’s law which takes the formwhere is the volumetric flow rate through the pore throat connecting pore bodies and , is the length of the pore throat, is the cross-sectional area, is the viscosity, and and are the pressures in the corresponding pore bodies. The above relationship is used for every pore throat to obtain the flow rates and therefore the absolute permeability may be obtained. For two-phase flows, things become more complicated by the existence of the interfaces and the capillary effects (El-Amin et al. [78], Naraghi and Javadpour [79], Zhang et al. [80, 81], Lia et al. [82], Landry et al. [83], Gerami et al. [84], and Dong et al. [85]). A good review can be found in the work of Joekar-Niasar and Hassanizadeh [86].

9. Macroscopic Governing Equations

The permeability of shale reservoirs is very low compared to the permeability of conventional reservoirs. The gas production from a shale reservoir depends on the existence of natural or artificial fracture networks. So, the mathematical models describing flow and transport in fractured porous media are considered as the main framework to model fluid transfer between shale matrix and fractures. Traditionally, dual-continua models have been used to describe the transport in fractured porous media which consist of matrix blocks and fractures. In these models, Darcy’s law is assumed to be the main flow driver. The transport in matrix blocks has been described by several mechanisms such as Knudsen diffusion in nanopores, desorption from kerogen, and diffusion in solid. For example, Hadjiconstantinou [87] and Javadpour [46] discussed different mechanisms to investigate the shale-gas production rates. They employed Klinkenberg effect which considers the slippage coefficient to describe the diffusion mechanism. On the other hand, during depressurization of shale-gas reservoirs, induced stresses on the rock system can further reduce the pore space, which affects permeability. To account for such stress-dependent permeability upon depressurizing the reservoir, Raghavan and Chin [88] proposed a formula that considers permeability as a function of the pressure. Other researchers such as Raghavan and Chin [88] and Chipperfield et al. [89] proposed other relationships defining how permeability is related to the stress state. Other modeling approaches including fractional derivative formulation have also been implemented (e.g., [90]). The modeling of shale-gas reservoir is divided into two different main models. The first one is the dual-continua model (e.g., [80, 81, 91]) and the second one is the discrete fracture model (e.g., [92, 93]).

10. Dual-Continuum Model

Warren and Root [94] were the first to develop an idealized model to study the behavior of flow in fractured porous media based on dual-continua models. Bustin et al. [31] used a standard dual-continuum model to study the permeability effect. Ozkan et al. [95] used the dual-porosity model, in which the matrix was considered as uniform radius spherical blocks. They considered both the matrix diffusive Darcy flow and the fractures stress-dependent permeability for naturally fractured reservoirs ignoring sorption and desorption processes. Also, Moridis et al. [96] used the standard dual-continuum model to describe several mechanisms in kerogen. Wu and Fakcharoenphol [97] used a generalized dual-continuum methodology and proposed general reservoir simulators ignoring the adsorption and desorption processes. They implemented a general theoretical fracture model for simulating fluid and heat flow processes in fractured unconventional reservoirs. Guo et al. [98] developed a mathematical model including the above-mentioned mechanisms to describe the flow behavior in tight shale-gas formation. Recently, El Amin [99] presented an analytic solution using the power-series method for the apparent permeability Klinkenberg model. Arbogast and his coworkers made significant contributions in deriving the dual-porosity models (e.g., [100103]). Furthermore, Showalter and his coworker highlighted some of the mathematical properties of the dual-porosity models (e.g., [104106]).

In the matrix blocks of shale strata, free gas and absorbed gas coexist with each other. Mass accumulation term for free gas per unit volume is , for a single-phase gas reservoir. The adsorbed gas is estimated as 20%85% of shale gas. Mass accumulation term which describes adsorbed gas on the matrix surface is [107], where is the adsorbed gas volume per unit area of shale surface. The most common way to describe this process is the Langmuir isotherm model [18, 29, 34, 40], which is expressed as where is the mole volume under standard condition ; is the Langmuir volume; is the Langmuir pressure; is the density of shale core; is the matrix pressure; is the molecular weight. The two mass accumulation terms of both free and absorbed gas are combined, . The real gas law is considered, , where p is the pressure, is the gas constant, is the temperature, is the mass, is the gas volume, and is the gas deviation factor. The mass density of the gas becomes . In order to calculate , one may use the cubic Peng-Robinson equation of state [108]:where , , , and . and are the critical temperature and critical pressure, respectively. The Klinkenberg method [20] was used to correct the effective gas permeability to a liquid-equivalent permeability using a “gas-slippage” factor. The Klinkenberg effect becomes significant in modeling gas flow in reservoirs in systems with low pressure or low permeability. The matrix apparent permeability can be written as [46] where , where is the tangential momentum accommodation coefficient, which takes values within the range . For the gas flow in fractures, we consider Knudsen diffusion and viscous flow. The mass flux may be represented as where . The apparent permeability of fractures is expressed as Knudsen diffusion coefficient, , for the fracture system is defined as [46] where is the initial fracture porosity. The dual-porosity dual permeability (DPDP) model has two mass conservation equations, one for the matrix blocks continuum and the other for the fractures continuum, coupled with a transfer term. The mass exchange between matrix blocks and fractures is represented by a shape factor [109]. Warren and Root [94] defined the shape factor for cubic matrix blocks as , where is the set of normal fractures, is a characteristic length, and , and are lengths of the sides of a cubic matrix block. Some other methods have been considered to handle the matrix-fracture connection such as the boundary conditions method [110]. The transfer of gas between the matrix and fracture systems is represented by where is the crossflow coefficient between the fracture and matrix systems. The fluid flows into the fracture from the matrix represented by a sink term. In order to define this sink term, we present the model developed by Aronofsky and Jenkins [111] for gas production from a vertical well; namely, where when the production well is placed in the center, while when the production well is located in the corner; is the bottom hole pressure; is the average fractures pressure around the well; is well radius; is the drainage radius, which can be calculated by where and are the grid lengths in and directions. Terzhagi (1936) reported that rock deformation based on effective stresses has also some effects on the transport of fluids in fractured systems. Biot and Willis [112] developed a generalization for the effective stress model such that the permeability and the porosity are both dependent on the stress. The porosities are related to the mean effective stress according to the following set of equations [97]: where is the effective stress which is defined by ; is the mean total stress in matrix blocks; is the mean total stress in fractures; are Biot’s effective parameters. The derivative of the porosity is given as where is the initial porosity of the matrix system; is the initial porosity of the fractures system. The stress evaluation is added to each term of the mass accumulation equation. Therefore, the DPDP model consists of two equations that can be finally written aswhere and .

11. Discrete Fracture Model

The discrete fracture model (DFM), a well-known reduced modeling technique for the simulation of flow and transport in the fractured porous media [113115], is reviewed in this section. In the single discrete fracture modeling, each fracture is represented explicitly using high resolution mesh (El Amin et al. [116]). So, as fractures have high permeability and big variations in physics, a denser mesh is required in fractures than in the matrix. For example, Baca et al. [117] presented a two-dimensional model for single-phase flow with heat and solute transfer in fractured porous media. Juanes et al. [118] proposed a finite-element formulation for a single-phase flow in fractured formations. Matthai et al. [119] introduced numerical investigations for two-phase flow in fractured porous media using control-volume finite-element method. Cipolla et al. [120] introduced numerical modeling using automated unstructured gridding method to simulate the well performance from the complex fractures. Sheng et al. [121] employed the extended finite-element method to study the gas transport of shale in a complex fracture network. The discrete fracture model for a single-phase flow with fluid exchange between the fracture and the surrounding rock matrix where fractures are treated as interfaces of dimension () was considered by Martin et al. [122]. Jaffré et al. [123] extended this model for the case of two-phase flow. The interaction between the matrix and the fracture is effected by Robin boundary condition along the two sides of the fracture, and the flux discontinuity through the fracture is represented by a source term. Also, nonlinear transmission conditions may be of relevance when considering reactive flow in fractured media. Pop et al. [124] employed nonlinear transmission conditions for reactive flow in fractured media. Recently, the phase field modeling of flow in fractured media has been considered by Lee et al. [125] and Mikelić et al. [126, 127]. The embedded discrete fracture models have been employed by Li and Lee [128] and Moinfar et al. [129] discretized the complex fractures into a number of segments while the matrix is treated as structured grids. Jiang and Younis [130] used the two hybrid methods, namely, embedded discrete fracture models with multiple interacting continua and the coupling of unstructured DFM with continuum type. The fracture networks’ complexity has important effects on the well performance. The fractures may be planar or nonplanar, natural or hydraulic. Some examples of different complex planar/nonplanar hydraulic/natural fracture networks are shown in Figure 11 [92]. The hydraulic fractures are connected with natural fractures and wellbores are connected in different schemes.

As the fracture network consists of nodes connected with segments, the pressure drop from the node node to the node node may be given as [131] where is the gas flow rate at node of the fracture segment, is the coefficient of Darcy flow for gas at the fracture segment , and is the coefficient of non-Darcy flow. is the gas flux at the fracture segment j, is the fracture permeability, is the fracture width, is the fracture height, and is the non-Darcy Forchheimer coefficient. Considering time variation from the time step to , the mass balance in the pore is given as the summation of all the mass flux connected with the pore [80, 81]:where is the slippage coefficient and is the pore center distance between the pore and the pore .

12. Numerical Methods for Solving the Governing Laws

On the other hand, simulation of transport phenomena is computationally challenging because it demands high accuracy and local mass conservation. Transport in geological media involves long durations; even tiny errors in each step can accumulate to a huge error. Porous media manifest dramatically different spatial and temporal scales. Heterogeneity, anisotropy, and discontinuity of medium properties require special treatment for computationally efficient approximation of advection, diffusion, dispersion, and chemical reactions. Some numerical approximation schemes fail to preserve important physical and/or mathematical principles and lead to erroneous simulation results. These challenges are addressed, for example, by Moortgat et al. [132], Sun and Wheeler [133, 134], Kou and Sun [135], Dawson et al. [136], Kou and Sun [137], Radu et al. [138], Radu and Pop [139], Vohralík and Wheeler [140], Vohralík [141], Marchand et al. [142], and many others. In general, using mixed finite-element method (MFEM) and finite volume methods in reservoir simulation and transport in porous media is required as they are locally mass conservative. On the other hand, new algorithms have been developed to speed up the computations that could, otherwise, be considerably slow. In particular, in programming languages that require repeating interpretations (like Matlab and Python), loops and logical statements can consume a considerable amount of time just for translation. Sun et al. [143] developed a technique that vectorizes all the difference operations replacing totally the loops. This makes the algorithm’s speed compare very well with those developed using languages that do not require repeated interpretation (like FORTRAN, C, and C++). Later, Salama et al. [144] generalized this algorithm to the problem of two-phase flows in porous media. Salama et al. [145, 146] also developed what they called the experimenting field algorithm in which the matrix of coefficients of a linear system is constructed automatically, rather than being inputted manually. This technique has been tested in various cases and shows to be very effective (Salama et al. [147149], Wu et al. [150], Negara et al. [151, 152], and El-Amin et al. [153]). Furthermore, Salama et al. [154] and Salama [155] developed a technique that numerically solves the problem of flow towards a well in an efficient manner. This technique is very accurate even when the mesh closer to the well is coarse.

The next paragraphs provide only a fleeting study of an otherwise rich and active area of research for numerical methods for flow and transport in porous media. For more discussions on this subject, the readers are advised to refer to excellent textbooks on this subject (e.g., Helmig [156], Chen et al. [157], and Nordbotten and Celia [158]). The conservation laws that describe various transport phenomena in porous media (e.g., those presented in this paper) are in the form of partial differential equations (PDEs). The analytical solutions of these equations are very limited to simple problems, geometries, and boundary conditions. Therefore, it is likely that researchers will look for numerical methods to solve various transport problems in porous media. In these methods, the solutions are obtained at discrete points inside the domain rather than at every point as is the case in analytical methods. Several classes of numerical methods have been adapted to solving different problems in porous media including finite differences, finite elements, finite volume, boundary elements, and spectral methods. The methodologies and requirements for each method can be classified into two general categories, namely, methods that directly approximate the derivatives (e.g., finite differences, classes of finite volume) and methods that approximate the solution functions themselves. In the former methods, the derivatives must exist to at least as much as the governing differential equations require. These methods, therefore, do not allow possible approximate solutions that do not satisfy the governing PDEs strongly. The later class of numerical methods have the advantage that they can generate approximate solutions that can satisfy the governing equations weakly, or in other words in an integral sense. The first methods are straightforward and therefore in this review we highlight some of the features of the second class of numerical methods.

The mixed method is based on the conservation law together with the constitutive law for the flux in terms of pressure (Darcy’s law) and solving the flux variable and the pressure together. Since this preserves the structure of the equations, this approach is in consistence with the required properties of the numerical methods. Though, in the previous section, equations are derived to model the flow using non-Darcy models, we limit the discussion here to the Darcy model. The mixed method being discussed here can be directly adapted to those developed in the previous section for nonzero . The weak formulation for the above equations is obtained by using , v as test functions, multiplying the governing equations of mass conservation and Darcy’s law, respectively, and integrating over the domain . This yieldsfor all , v. In the above equations, the notation implies the inner product which is defined asNote that, in the above, partial integration has been used for putting the derivative from pressure to test function v. Further, it has been assumed that the boundary conditions are such that the boundary terms cancel out. The inner product above is a standard inner product and the integration has been performed over the domain Ω. The test functions and v belong to different spaces justifying the nomenclature of the mixed method.

The finite-element spaces for the solutions and and for the test functions v and are chosen to be the same. However, the choice for vector spaces for and cannot be made independently and has to satisfy a compatibility condition known as the inf-sup condition. In general, the heterogeneities and the nonlinearities including advancing fronts in the porous media imply that the regularity of the solution is quite low. This motivates the use of low order finite-element space. Raviart and Thomas [159], later extended by Nédélec [160], introduced the first family of mixed finite-element spaces for the pressure equation under consideration. For the lowest order, approximation satisfying an inf-sup condition is given in terms of Raviart-Thomas elements, denoted as RT0. For this approximation, the scalar spaces (e.g., for pressure, ) are chosen as constants in each element of the mesh of the domain, whereas the vector spaces, say for the flux, , are chosen to be piecewise linear. The degrees of freedom for the flux unknowns are at the mesh edges. The flux unknowns are taken to be constant at each edge and the linear extrapolation in the interior of each element is taken to ensure that This ensures the inf-sup compatibility condition.

In the coupled flow and transport problems being considered here, the mixed method provides several advantages over the conformal finite-element method. It is noted that the heterogeneities in the permeability imply that the gradient of the pressure will have discontinuities and will be rough. However, the flux is smoother. Thus, the mixed method computes the flux explicitly and accurately and the flux has a physical meaning. In fact, the discrete equation for the flux represents the conservation equation and hence reflects the physics of the problem. In the transport model, the flux is used for the convection of the reactant species and the accuracy of the flux is an important consideration for the accuracy of concentration profiles. The immediate disadvantages are the increase in the number of unknowns and the linear system losing its positive definiteness. To overcome these deficiencies, a quadrature is used to approximate the flux unknowns in terms of pressure unknowns that can be locally inverted. The flux unknowns are then eliminated giving rise to a system with only pressure unknowns. Such a system is positive definite and it can be shown theoretically that the approximation in using the quadrature does not lead to any loss in convergence order [161].

13. Characteristic Method for Transport Equation

The multiphase flow or the transport of chemical species is characterized by the presence of front of the chemical species being transported or of saturation. This is due to the hyperbolic character of the governing equations and is manifested because of stronger advection compared to diffusion. In the multiphase flow, the behavior is described by the degeneracy of the governing models where the nature of equation changes from parabolic to hyperbolic, and in the transport of the chemical species case, the stronger advection including limited diffusion means that the numerical schemes have to include both features of parabolic and hyperbolic. However, the schemes used for parabolic equations do not work well when the governing equations are advection dominated, especially capturing the sharp fronts while ensuring the conservation of local mass. The design of numerical schemes therefore has to incorporate the hyperbolic nature of the transport equation for discretization. Several numerical methods have been developed that provide the adaptation from parabolic to hyperbolic: explicit method of characteristics, upstream-weighted finite differences [162], interior penalty Galerkin schemes [163], higher order Godunov schemes [164], streamline diffusion method [165], and modified method of characteristics-Galerkin finite-element procedure [166, 167]. The reader is referred to the textbook of Chen [168] for further discussions on the use of some of these techniques in porous media context.

14. Deformable Porous Media

The theoretical basis for the study of flow in deformable porous media based on the method of volume averaging can be found in the early work of Whitaker [169]. An important feature in tight and shale formations is that the property of the porous medium may evolve as a result of chemical reactions. Along with porosity, these properties include transport properties of the medium such as permeability and diffusivity. A nonnegligible change in the porosity may happen due to chemical reactions such as precipitation and dissolution [170]. The transport properties are typically considered as a function of porosity ([171175] and many others) but in an evolving pore skeleton, these properties may change considerably as they depend on the details of the pore-scale geometry. To consider the effects due to the evolving microstructure is of crucial importance for obtaining reliable upscaled models, because otherwise features like pore clogging or damaging of the structure will not be captured at all. The existing models for reactive flow either do not include these features or are restricted to very simple geometries (e.g., 1D) as in the work of Alshawabkeh and Rahbar [176], or they use ad hoc proposed laws based on the porosity-permeability relationship. Nevertheless, to include the evolution of the microscale is complicated due to the complexity of the domain at the pore scale and the occurrence of free and moving boundaries (pores have a variable, solution dependent structure, which is not known a priori). Recently, problems with an evolving microstructure were considered for simplified geometries and saturated flow (see [177] and van Noorden et al., 2010 [178, 179]). Extensions to multiphase models having an evolving microstructure coupled to saturated/unsaturated flow (i.e., the pore space may be filled with water and air) are quite open. A consistent theory for the derivation of deformable porous media taking into account the detailed pore-scale geometry evolution and the corresponding changes in the transport properties is by far incomplete. The macroscale models consist of fully coupled, nonlinear, degenerate parabolic partial differential equations for flow, reactive transport and heat, and additional ordinary differential equations for the microstructural evolution. Two degeneracy types appear here: parabolic/elliptic due to saturated/unsaturated flow and parabolic/hyperbolic due to pore clogging. The design and implementation of efficient and reliable numerical schemes for the macroscale system is therefore a very demanding task. The evolving geometry is reflected in the upscaled models through additional equations describing the topological changes describing the porosity and permeability. How to properly account for these topological changes is at present open despite early attempts made in the references cited above. For such models, it remains open to develop, implement, and analyze efficient multiscale, mass-conservative solution schemes. An immediate difficulty is in the increase of dimensions of the problem due to the fine scale evolving geometry description coupled to Darcy scale models.

Upscaling methods have thus far been mainly applied to rigid porous media. Any approach for upscaling should include the variation of the geometry at the pore scale. Recently, an extension of homogenization techniques was introduced by van Noorden [177] including a level set approach to capture the moving interfaces at pore scale. Assuming the fixed geometry for the precipitation-dissolution models allows performing rigorous analysis of the mathematical models. The related effective models describing convective-diffusive transport for the periodic case have been derived by Kumar et al. (2016). Effective equations for the convection dominated case leading to the Taylor dispersion including the moving boundary were developed by Kumar et al. [178, 180184].

There is no doubt that there are potential environmental risks associated with acquiring energy resources from unconventional reservoirs using fracking technologies. The societal dilemma, however, lies in how to devise methodologies for obtaining the natural gas while limiting environmental damage. Global warming, whose adverse effects are, nowadays, signified and felt throughout the world, represents a pressing issue. Escaping methane, probably, represents the most serious and noticeable impact because of its direct and acute effect on the local environment. There have been several instances where groundwater and even surface water get contaminated with the release of gases from shale formations in relatively larger quantities. Groundwater distributed to houses in the nearby areas of drilling sites has been shown, in some cases, to catch fire due to the release of gases [185187]. Such release of gases (mainly methane) not only contaminates water bodies but also contributes to air pollution. Global warming potential (GWP), a number that allows comparing the potential impact of different greenhouse gases (GHG) on global warming, indicates that methane warms the planet 86 times as much as CO2, according to the Intergovernmental Panel on Climate Change. Methane is 25 times more potent in trapping heat in the atmosphere than carbon dioxide. The major source of gas release from unconventional reservoirs may, probably, be attributed to the fracking technology used. That is, there is no guarantee that the stimulated fractures will be confined to where the wellbore exists. In other words, there are no measures to ensure that the stimulated fractures do not extend to the top edge of the formations. If this happens, the gases can probably find a pathway to the top layer of the formation, thereby to the local groundwater reservoirs, and finally to the atmosphere. The risks and concerns of hydraulic fracking are generally associated with the contamination of groundwater, methane pollution and its impact on climate change, air pollution impacts, exposure to toxic chemicals, blowouts due to gas explosion, waste disposal, large volume water use in water-deficient regions, fracking-induced earthquakes, workplace safety, and others. In addition, chemical additives that are used in the drilling mud, slurries, and fluids, which are essentially required for the fracking processes, can generally leak into the local environments through fissures, sealing, cracks, faulty design, or construction of the cement well casings, among others. Each well utilizes millions of gallons of toxic fluid containing not only the added chemicals, but also other naturally occurring radioactive materials, liquid hydrocarbons, brine water, and heavy metals (Osborn et al. [188], Bamberger and Oswald [189], and Colborn et al. [190]). Fissures created by the fracking process can also create underground pathways for gases, chemicals, and radioactive materials.

16. Discussion and Conclusion

This paper reviewed aspects related to the flow and transport phenomena in tight and shale formations. These formations are believed to contain a considerable amount of hydrocarbon resources. The major characteristics of such formations compared with conventional reservoirs are that their permeability is so small that traditional techniques for producing oil and gas may not work. Unconventional methods are, therefore, required to be developed to facilitate the transport of oil and/or gas stored in these formations. Horizontal wells have been shown to be very effective in increasing the area of contact between the wellbore and the formation. To facilitate the transport of hydrocarbons, fractures are formed using fracking technology. Such fractures expose a larger area of the rock mass to thermodynamic conditions which results in the mobilization of entrapped and adsorbed gas. There are several aspects related to the process of fracturing the rock mass. These are related to the fact that there is no guarantee that these induced fractures are confined to the vicinity of the well leg. In other words, fractures may extend to the top layer of the formation, providing the released gas with a pathway to the local environment. There are several aspects related to the mechanisms involved in the release and transport of stored gas towards the production wells. Micro-CT analysis of rock samples indicates that the gas exists mostly in an organic region of the formation called kerogen. The gas is found in these regions adsorbed to the surface. When the pressure is decreased, the adsorbed molecules are derived to the pore space and then to the fractures network. The sizes of pore space in such formation give rise to whether the continuum hypothesis may be applicable in such small size conduits. Indeed, the characteristic length scale of pore space is approximately in the same order of magnitude as that of the average distance between gas particles. The Knudsen number, which compares the mean free path of the molecules and the characteristic length scale of the flow conduits, was shown to be relatively large, violating therefore the constraints set by the continuum hypothesis. This necessitates the need to consider multiscale approaches including molecular scale, pore scale, and continuum scale phenomena. The fact that the Knudsen number may be larger than one highlighted the need to consider interesting physical phenomena, including flow slippage and Knudsen diffusion. Under these scenarios, the permeability has been shown to be a function of the pressure. Furthermore, the depressurization of shale-gas reservoirs induces several stress-induced deformations that, likewise, affect the permeability. To model these complex processes at the macroscopic level, two general frameworks have been proposed: the dual-continua approaches and the discrete fractures model. Several numerical techniques have been used in solving the governing equations including mixed finite-element methods, finite volume methods, and the method of characteristics. A comprehensive list of recent references have been included to help researchers interested in pursuing research in this field to find the appropriate material.

Conflicts of Interest

The authors declare that they have no conflicts of interest.


Shuyu Sun acknowledges that this work is supported by the KAUST research fund awarded to the Computational Transport Phenomena Laboratory at KAUST through Grant BAS/1/1351-01-01. In addition, Shuyu Sun thanks his Ph.D. student Tao Zhang for his help in collecting and compiling results on Klinkenberg effect and Knudsen diffusion in shale.