Review of Formation and Gas Characteristics in Shale Gas Reservoirs

: An accurate understanding of formation and gas properties is crucial to the e ﬃ cient development of shale gas resources. As one kind of unconventional energy, shale gas shows signiﬁcant di ﬀ erences from conventional energy ones in terms of gas accumulation processes, pore structure characteristics, gas storage forms, physical parameters, and reservoir production modes. Traditional experimental techniques could not satisfy the need to capture the microscopic characteristics of pores and throats in shale plays. In this review, the uniqueness of shale gas reservoirs is elaborated from the perspective of: (1) geological and pore structural characteristics, (2) adsorption / desorption laws, and (3) di ﬀ erences in properties between the adsorbed gas and free gas. As to the ﬁrst aspect, the mineral composition and organic geochemical characteristics of shale samples from the Longmaxi Formation, Sichuan Basin, China were measured and analyzed based on the experimental results. Principles of di ﬀ erent methods to test pore size distribution in shale formations are introduced, after which the results of pore size distribution of samples from the Longmaxi shale are given. Based on the geological understanding of shale formations, three di ﬀ erent types of shale gas and respective modeling methods are reviewed. Afterwards, the conventional adsorption models, Gibbs excess adsorption behaviors, and supercritical adsorption characteristics, as well as their applicability to engineering problems, are introduced. Finally, six methods of calculating virtual saturated vapor pressure, seven methods of giving adsorbed gas density, and 12 methods of calculating gas viscosity in di ﬀ erent pressure and temperature conditions are collected and compared, with the recommended methods given after a comparison.


Introduction
Shale gas has received great attention from governments all over the world, especially after the successful shale gas revolution in North America. As a type of clean energy with abundant reserves, shale gas is believed to be one of the most promising replacements for conventional energy in the future. According to the Energy Information Administration (EIA) (2016), shale gas production is expected to drive world natural gas production growth in the coming decades and will account for approximately 30% of world natural gas production by 2040 (Figure 1a). The United States and China are predicted to be the two largest shale gas producers in the world by the end of the forecast period, with shale gas production making up 70% and 40%, respectively, of each country's total natural gas production (Figure 1b). Taking China as an example, due to the relatively high economic growth and increasing attention to environmental protection, natural gas consumption is expected to increase from 19 Bcf/d in 2015 to 57 Bcf/d in 2040 (Figure 1c), accounting for a quarter of all global natural gas consumption growth between 2015 and 2040 (EIA, 2017). Driven by the development of shale gas resources, China's domestic natural gas supply will grow from 13 Bcf/d in 2016 to 39 Bcf/d by 2040, with shale gas production increasing from 0.7 Bcf/d in 2016 to 10 Bcf/d by 2030 and 19 Bcf/d by 2040. As we can see from Figure 1d, shale gas is expected to increase the fastest and will account for more than 30% of the total natural gas supply in China by 2040.
As one of the unconventional energy resources, shale gas reservoirs have uniqueness and complexity in terms of the gas storage type, transporting mechanism, and reservoir development mode, which makes the commercial production of shale plays a very challenging task for petroleum engineers. The United States, Canada, and China are the only three countries that produced commercial volumes of natural gas from shale formations by 2015. Although horizontal well drilling and hydraulic fracturing have been applied to produce shale gas in Australia and Russia, no commercial gas volumes were obtained from low-permeability shale formations. Currently, the commercial development of shale gas resources in North America and China mainly benefits from advanced engineering technology. The theoretical understanding of shale gas storage capacity, gas transporting mechanism in nanopores or micropores, and pore structure characterization is still not clear, being far behind engineering practice [1].
In this review, our attention will mainly be paid to three aspects: (1) petrological, organic geochemical characteristics and micropore structures of shale formations; (2) different types of adsorption models as well as their principles and application range, including Gibbs excess sorption, supercritical adsorption phenomenon, and adsorption/absorption models; (3) different methods of calculating gas physical properties, such as virtual saturated vapor pressure, adsorbed gas density, free gas density, free gas viscosity, etc. Different models on each subject will be compared and evaluated based on their physical meaning, reliability, accuracy, and applicability, which are significant for accurate numerical simulation and enhancing hydrocarbon recovery in shale gas reservoirs.  Figure 1d, shale gas is expected to increase the fastest and will account for more than 30% of the total natural gas supply in China by 2040.
As one of the unconventional energy resources, shale gas reservoirs have uniqueness and complexity in terms of the gas storage type, transporting mechanism, and reservoir development mode, which makes the commercial production of shale plays a very challenging task for petroleum engineers. The United States, Canada, and China are the only three countries that produced commercial volumes of natural gas from shale formations by 2015. Although horizontal well drilling and hydraulic fracturing have been applied to produce shale gas in Australia and Russia, no commercial gas volumes were obtained from low-permeability shale formations. Currently, the commercial development of shale gas resources in North America and China mainly benefits from advanced engineering technology. The theoretical understanding of shale gas storage capacity, gas transporting mechanism in nanopores or micropores, and pore structure characterization is still not clear, being far behind engineering practice [1].
In this review, our attention will mainly be paid to three aspects: (1) petrological, organic geochemical characteristics and micropore structures of shale formations; (2) different types of adsorption models as well as their principles and application range, including Gibbs excess sorption, supercritical adsorption phenomenon, and adsorption/absorption models; (3) different methods of calculating gas physical properties, such as virtual saturated vapor pressure, adsorbed gas density, free gas density, free gas viscosity, etc. Different models on each subject will be compared and evaluated based on their physical meaning, reliability, accuracy, and applicability, which are significant for accurate numerical simulation and enhancing hydrocarbon recovery in shale gas reservoirs.  [2][3][4].
In this section, we summarize the petrological and organic-geochemical characteristics of shale gas reservoirs in Sichuan basin in China compared to the shale gas development in North America. The chosen samples are marine shale of the Lower Silurian Longmaxi Formation. The porosity and permeability of the main shale gas reservoirs in North America and China are collected and tested, respectively, based on which the storage space types and pore size distributions are analyzed. Finally, different kinds of shale gas as well as their modeling methods are identified and compared.

Petrological and Geochemical Characteristics
In this part, our previous work related to shale gas reservoir characterization and assessment is introduced. Samples of the Longmaxi Formation in south Sichuan basin are collected as experimental objects.
The highly mature Longmaxi marine shale is one of the most important candidates for the commercial development of shale gas resources in China [10]. The total organic carbon (TOC) content ranges from 0.4% to 18.4%, with the organic matter (OM) mainly composed of type I and II1 kerogen [10,11]. Its vitrinite reflectance (R0) values range from 1.8% to 4.2% [11]. The Longmaxi Shale is found to be porous and permeable [11,12], with porosity ranging from 1.2% to 10.8% and permeability ranging from 0.25 µD to 1.737 mD. Other geological and petrophysical characteristics of the Longmaxi Shale Formation can be found in previous publications [13,14].  [2][3][4].
In this section, we summarize the petrological and organic-geochemical characteristics of shale gas reservoirs in Sichuan basin in China compared to the shale gas development in North America. The chosen samples are marine shale of the Lower Silurian Longmaxi Formation. The porosity and permeability of the main shale gas reservoirs in North America and China are collected and tested, respectively, based on which the storage space types and pore size distributions are analyzed. Finally, different kinds of shale gas as well as their modeling methods are identified and compared.

Petrological and Geochemical Characteristics
In this part, our previous work related to shale gas reservoir characterization and assessment is introduced. Samples of the Longmaxi Formation in south Sichuan basin are collected as experimental objects.
The highly mature Longmaxi marine shale is one of the most important candidates for the commercial development of shale gas resources in China [10]. The total organic carbon (TOC) content ranges from 0.4% to 18.4%, with the organic matter (OM) mainly composed of type I and II1 kerogen [10,11]. Its vitrinite reflectance (R 0 ) values range from 1.8% to 4.2% [11]. The Longmaxi Shale is found to be porous and permeable [11,12], with porosity ranging from 1.2% to 10.8% and permeability ranging from 0.25 µD to 1.737 mD. Other geological and petrophysical characteristics of the Longmaxi Shale Formation can be found in previous publications [13,14].

Mineral Composition
The rock mechanics, adsorption capacity, and well productivity of shale gas reservoirs are directly determined by the relative content of different minerals due to the property differences among the minerals. Taking 18 samples from X1 well, eight samples from X2 well, and 60 samples from X3 well, we tested the mineral compositions using an X'Pert Pro type X-ray scattering diffractometer produced by PANalytical B.V. (Almelo, Netherlands ) and following the Chinese Oil and Gas Industry Standard SY/T 5163 1995 and SY/T 5983 94. The laboratory temperature is 24 °C and the humidity is 30%. The test results are shown in Figure 3. As shown, the main mineral compositions of the Longmaxi Formation are quartz (11-70%, average 31.06%) and clay minerals (7-64%, average 33.87%). Comparing the results with those in the Mississippian Barnett Shale of the Fort Worth Basin in North America [15,16], the content of brittle minerals (such as quartz, feldspar, and calcite) in the Longmaxi Formation is relatively lower, while the content of clay minerals is higher (Figure 4).   [15,16].  The rock mechanics, adsorption capacity, and well productivity of shale gas reservoirs are directly determined by the relative content of different minerals due to the property differences among the minerals. Taking 18 samples from X1 well, eight samples from X2 well, and 60 samples from X3 well, we tested the mineral compositions using an X'Pert Pro type X-ray scattering diffractometer produced by PANalytical B.V. (Almelo, Netherlands ) and following the Chinese Oil and Gas Industry Standard SY/T 5163 1995 and SY/T 5983 94. The laboratory temperature is 24 °C and the humidity is 30%. The test results are shown in Figure 3. As shown, the main mineral compositions of the Longmaxi Formation are quartz (11-70%, average 31.06%) and clay minerals (7-64%, average 33.87%). Comparing the results with those in the Mississippian Barnett Shale of the Fort Worth Basin in North America [15,16], the content of brittle minerals (such as quartz, feldspar, and calcite) in the Longmaxi Formation is relatively lower, while the content of clay minerals is higher (Figure 4).     With less clay mineral and more brittle minerals, natural or induced fractures more easily develop under external forces. On the contrary, the higher clay mineral content has a negative effect on volume stimulation, since most of the energy is absorbed by shale formations. Under such circumstances, Energies 2020, 13, 5427 5 of 50 plane fractures are more likely to be generated, rather than tree-like or reticular structural fractures. Generally, the brittle minerals need to be higher than 40% and the clay minerals need to be less than 30% for a potential shale gas reservoir to be commercially developed [16,17].
The content of different minerals also exhibits different trends in terms of formation depth, as shown in Figure 5. For the formation depth increasing from 2120 m to 2250 m, the content of calcite decreases from 12% to 5%, and the content of brittle minerals increases from 28% to 60%. The content of clay minerals increases from 44% to 48% for the formation depth, increasing from 2120 m to 2205 m, while it decreases rapidly from 48% to 8% for formation depth increasing from 2205 m to 2250 m.
Energies 2020, 13, 5427 5 of 51 With less clay mineral and more brittle minerals, natural or induced fractures more easily develop under external forces. On the contrary, the higher clay mineral content has a negative effect on volume stimulation, since most of the energy is absorbed by shale formations. Under such circumstances, plane fractures are more likely to be generated, rather than tree-like or reticular structural fractures. Generally, the brittle minerals need to be higher than 40% and the clay minerals need to be less than 30% for a potential shale gas reservoir to be commercially developed [16,17].
The content of different minerals also exhibits different trends in terms of formation depth, as shown in Figure 5. For the formation depth increasing from 2120 m to 2250 m, the content of calcite decreases from 12% to 5%, and the content of brittle minerals increases from 28% to 60%. The content of clay minerals increases from 44% to 48% for the formation depth, increasing from 2120 m to 2205 m, while it decreases rapidly from 48% to 8% for formation depth increasing from 2205 m to 2250 m. The TOC content is determined by a CS230 carbon/sulfur analyzer (LECO, St. Joseph, USA) with samples crushed into powder less than 100-mesh. The powder is pyrolyzed up to 600 °C and the inorganic carbon is removed by hydrochloric acid. The relationship between quartz and TOC content can help to explain the origin of quartz, i.e., detrital or biogenic.
There is no linear relationship between detrital quartz and TOC content, while a positive correlation can be found for biogenic quartz and TOC. As shown in Figure 6, the correlation coefficients between quartz and TOC content in X1, X2, and X3 wells are larger than 0.53, indicating that the quartz in the targeted formation is biogenic. The TOC content is determined by a CS230 carbon/sulfur analyzer (LECO, St. Joseph, USA) with samples crushed into powder less than 100-mesh. The powder is pyrolyzed up to 600 • C and the inorganic carbon is removed by hydrochloric acid. The relationship between quartz and TOC content can help to explain the origin of quartz, i.e., detrital or biogenic.
There is no linear relationship between detrital quartz and TOC content, while a positive correlation can be found for biogenic quartz and TOC. As shown in Figure 6, the correlation coefficients between quartz and TOC content in X1, X2, and X3 wells are larger than 0.53, indicating that the quartz in the targeted formation is biogenic.

Organic Geochemical Characteristics of Shale
The organic matter richness, thermal maturity, and kerogen types are three key parameters for the accurate assessment of hydrocarbon-forming conditions. The organic matter richness not only affects  The organic matter richness, thermal maturity, and kerogen types are three key parameters for the accurate assessment of hydrocarbon-forming conditions. The organic matter richness not only affects the hydrocarbon generating strength, but also the development of organic pores and the adsorbed gas content. The lower limit value of the TOC content for economic exploitation of shale gas reservoirs is approximately at 2.5-3 wt % [18]. However, with the development of technology, this value could become even lower.
Measuring 122 samples of TOC content, we find that the TOC content ranges from 0.43% to 8.39%, with an average value of 2.20% in Longmaxi Formation of Sichuan Basin. The samples with a TOC content less than 2.00% account for 57.38%, while 42.62% of samples have a TOC content larger than 2.00%. This reflects the fact that TOC is abundant in Longmaxi Formation, which is advantageous for shale gas generation and storage. However, comparing with shale gas reservoirs in North America, the TOC content in Sichuan Basin is smaller. The TOC content of Antrim shale and New Albany shale is between 1% and 25%, while it is between 0.45% and 4.5% in the Barnett Shale and Lewis Shale [19].
Analyzing the TOC content data of three wells in the Longmaxi Formation longitudinally, we find a positive relationship between TOC content and depth, as shown in Figure 7. The TOC content at the bottom of the formation is much larger than that at the top. The TOC content increases with depth at 2.3-10.0% per 100 m in the targeted formation.
The kerogen types can be classified into sapropelic type (type I), mixed -type (type II), and humic type (type III) [20]. All three types of kerogen can generate natural gas. For type I and II1 kerogen, oil is generated first and then cracked into gas. The type III kerogen is not advantageous for oil generation and gas is formed directly from organic matter [17]. The abundance of organic matter is the material basis for hydrocarbon generation, while the type of organic matter determines the hydrocarbon generating potential and hydrocarbon characteristics.
The thermal evaluation extent of organic matter can be characterized by thermal maturity, reflected by the vitrinite reflectance R o , which is the basis of assessing the hydrocarbon generating potential of source rocks ( Table 1). The organic matter maturity range of 1.1% < R o < 3.5% is advantageous for the generation of shale gas [3,21]. Single well production in more mature shale formations is larger than in less mature ones, because more gas is generated by kerogen or the thermal cracking of crude oil in more mature areas. The average thermal maturity of the targeted Barnett Shale is R o = 1.7%, with the maximum value larger than 2.0% [22]. Oil and gas are generated from the kerogen with initial R o < 1.1% in the Barnett Shale, but gas is produced within the formation in the Newark East and surrounding areas at higher thermal maturity, i.e., R o > 1.1% [23]. The effect of organic matter maturity R o on shale gas reservoirs is very complicated and needs further study. TOC content in Sichuan Basin is smaller. The TOC content of Antrim shale and New Albany shale is between 1% and 25%, while it is between 0.45% and 4.5% in the Barnett Shale and Lewis Shale [19].
Analyzing the TOC content data of three wells in the Longmaxi Formation longitudinally, we find a positive relationship between TOC content and depth, as shown in Figure 7. The TOC content at the bottom of the formation is much larger than that at the top. The TOC content increases with depth at 2.3-10.0% per 100 m in the targeted formation. The kerogen types can be classified into sapropelic type (type I), mixed -type (type II), and humic type (type III) [20]. All three types of kerogen can generate natural gas. For type I and II1 kerogen, oil is generated first and then cracked into gas. The type III kerogen is not advantageous for oil generation and gas is formed directly from organic matter [17]. The abundance of organic matter is the material basis for hydrocarbon generation, while the type of organic matter determines the hydrocarbon generating potential and hydrocarbon characteristics.
The thermal evaluation extent of organic matter can be characterized by thermal maturity, reflected by the vitrinite reflectance Ro, which is the basis of assessing the hydrocarbon generating potential of source rocks ( Table 1). The organic matter maturity range of 1.1% < Ro < 3.5% is advantageous for the generation of shale gas [3,21]. Single well production in more mature shale formations is larger than in less mature ones, because more gas is generated by kerogen or the thermal cracking of crude oil in more mature areas. The average thermal maturity of the targeted Barnett Shale is Ro = 1.7%, with the maximum value larger than 2.0% [22]. Oil and gas are generated from the kerogen with initial Ro < 1.1% in the Barnett Shale, but gas is produced within the formation in the Newark East and surrounding areas at higher thermal maturity, i.e., Ro > 1.1% [23]. The effect of organic matter

Porosity and Permeability Characterization
Porosity and permeability are the two most important parameters to characterize gas storage and seepage capacities in shale gas reservoirs. Compared to conventional reservoirs, shale gas reservoirs are ultra-tight formations with extremely low porosity and permeability. Corresponding formation physical properties of main shale gas reservoirs in North America are attributed in Table 2 [3,15,18,19,25,26].  As the earliest country to commercially develop shale gas resources, the USA has formulated a standard system to evaluate the physical properties of shale formations. The Gas Research Institute (GRI) of America proposed a test method for shale cores to determine both the total porosity and the gas-bearing porosity of the shale matrix. Generally, the porosity range of shale formations is between 2% and 15%. From the statistical data in Table 2, we can see that the total porosity of shale formation in North America is between 2% and 14%, and the average value is between 4.22% and 6.51%. Following the procedure of GRI, the statistical results of measured gas-filled porosity is between 1% and 7.5%, and water-filled porosity is between 1% and 8% in the Longmaxi Formation. The permeability of measured samples in the Longmaxi Formation is less than 0.1 mD, and the average pore-throat radius is smaller than 0.002 µm.

Pore Structure Division
Pore structures in shale formations can roughly be divided into two types: matrix pores and fractures (Table 3). Matrix pores are the main storage space of shale gas, directly determining the reserve of a shale gas reservoirs. The development of fractures as well as the connectivity of pores determines the gas-transporting and -producing capabilities [27].

Major Classes Subgroups Types
Fractures Based on previous studies [16,[28][29][30][31], we came up with a new classification method of gas storage space in shale formations in this paper by comprehensively utilizing core observation, thin-section analysis, SEM analysis, and emission scanning electron microscopy after argon ion Energies 2020, 13, 5427 9 of 50 polishing technologies. The main storage and seepage space in shale formations can be classified into fractures and matrix pores (Table 3). According to the origin mode, matrix pores are further divided into inorganic pores and organic pores, among which inorganic pores include intergranular pores and intragranular pores. Fractures can be classified into tectonic extensional fractures, tectonic shear fractures, interlayer bedding fractures, rock convergent fractures, and abnormal pressure fractures (Table 3).

Pore Size Distribution and Influential Factors
Gas storage and seepage mechanisms vary significantly due to the difference of pore sizes. Clear knowledge of pore size distribution (PSD) in shale formations is essential for shale gas exploitation and development. MICP, gas adsorption, and NMR are three commonly used methods to determine PSD.

Pore Size Distribution and Influential Factors
Gas storage and seepage mechanisms vary significantly due to the difference of pore sizes. Clear knowledge of pore size distribution (PSD) in shale formations is essential for shale gas exploitation and development. MICP, gas adsorption, and NMR are three commonly used methods to determine PSD. In this section, the principles and results of different methods will be introduced and analyzed, based on the measurements of samples from the Longmaxi Formation in Sichuan Basin, China.

Test Methods and Principles N 2 Adsorption Measurement
The Brunauer-Emmett-Teller (BET) adsorption model [32] is adopted to determine the specific surface of shale samples when 0.05 < p/p s < 0.35. The pressure is too small to achieve multilayer adsorption when p/p s < 0.05, and capillary condensation may happen when p/p s > 0.40. Some studies (e.g., [33]) pointed out that capillary condensation could not happen in shale gas reservoirs, since the common shale gas reservoir temperature is much higher than the critical temperature of shale gas (mainly methane). The two-parameter BET equation can be expressed as follows [32]: where V is the adsorbed gas volume, mL; V m is the saturated adsorption volume of monolayer, mL; p is pressure, Pa; p s is saturated vapor pressure, Pa; and b is a dimensionless constant related to the adsorption capacity. After measuring the adsorbed gas amount G, a linear relationship between p/[V(p s − p)] and p/p s (0.05 < p/p s < 0.35) can be found. According to the slope and intercept of the straight line, the saturated adsorption amount V m can be calculated, by which the specific surface of samples can be obtained: where N A is the Avogadro constant; A m is the cross-section area of N 2 (0.162 nm 2 ); W is the weight of the samples, g; and S g is the specific surface of samples, m 2 /g. The Barrett-Joyner-Halenda (BJH) equation [34] is used to calculate PSD when p/p s > 0.40: where γ is the surface tension, N/m; R is the mole heat capacity, J/(mol·K); T is the environmental temperature, K; and r is the pore radius, m.

Mercury Injection Capillary Pressure
A mercury intrusion porosimeter is widely adopted to determine PSD in conventional sandstone reservoirs, where the pressure of mercury and pore radius r satisfy the Washburn equation [35]: where ξ is the contact angle between mercury and shale surface; σ is surface tension of mercury, 10 −3 N/m; and p is the injection pressure, Pa.
The smallest pore radius that can be tested is determined by the highest pressure that mercury porosimetry can hold. In our study, a PoreMaster 60 mercury porosimeter is employed and its measurement range of pore size lies between 3.6 nm and 950 µm, but values near the lower limit can hardly be detected. This is because it is difficult to inject mercury into micro/nanopores due to the Energies 2020, 13, 5427 11 of 50 high capillary pressure. Meanwhile, high pressure may create artificial crack and stress sensitivity, which reduces the credibility of the measurement. Therefore, MICP is mainly used to analyze mesopores and macropores in shale samples.

Nuclear Magnetic Resonance (NMR)
The relaxation characteristic of a hydrogen nucleus under an external magnetic field is used to obtain the PSD by the NMR method, which causes no harm to the shale samples. The speed of relaxation is characterized by the longitudinal relaxation time T 1 and transverse relaxation time T 2 . The relaxation characteristics of fluid in different-sized pores are different, based on which PSD can be calculated. The transverse relaxation time T 2 is composed of bulk phase relaxation T 2B , surface relaxation T 2S , and diffusion relaxation T 2D , which is expressed as follows: The diffusion relaxation speed can be ignored compared to the surface relaxation speed in a uniform magnetic field, and the reciprocal of diffusion relaxation time T 2D is almost 0. Meanwhile, the bulk phase relaxation time T 2B is much bigger than the surface relaxation time T 2S . Therefore, 1/T 2D and 1/T 2B in Equation (5) can be ignored, and we have Letting C = F s · ρ 2 , we obtain the relationship between relaxation time T 2 and pore radius r via Note that the transformation coefficient C in Equation (7) is an empirical parameter varying from one area to another, which can be determined by experiments [36]. In order to obtain the value of C, the T 2 spectrum need to be measured for specific shale samples first, and a N 2 adsorption test needs to be conducted on exactly the same samples (or samples from the same formation) afterwards [36]. Since core plugs are used in NMR and crushed rock samples are needed for the N 2 adsorption test, it is essential to perform the NMR test prior to N 2 adsorption. Comparing the two PSD results, the transformation coefficient C can be fitted. Accurate determination of C is a key part of determining the PSD of shale samples by NMR.

Pore Size Distribution N 2 Adsorption Results
According to the principle of N 2 adsorption on measuring PSD, the measurable pore-throat size range is on the magnitude of nanometers, mainly micropores (<2 nm) and mesopores (2-50 nm). N 2 adsorption results ( Figure 9) show that PSD displays a high single peak at the pore size range of 2 to 5 nm, implying that nanopores ranging from 2 nm to 5 nm are very developed in shale formations, which is advantageous for the storage of adsorbed gas. Meanwhile, we see that pores larger than 10 nm are not very developed. In a PSD frequency histogram ( Figure 10), mesopores account for the largest percentage of all pores, followed by micropores and macropores. Micropores and mesopores accounted for more than 90% of the total pore volume.
Energies 2020, 13, 5427 12 of 51 is advantageous for the storage of adsorbed gas. Meanwhile, we see that pores larger than 10 nm are not very developed. In a PSD frequency histogram ( Figure 10), mesopores account for the largest percentage of all pores, followed by micropores and macropores. Micropores and mesopores accounted for more than 90% of the total pore volume.

Mercury Intrusion Results
Two peaks can be found on the MICP measurement results of PSD, as shown in Figure 11. The left peak is relatively small and smooth, corresponding to macropores of 10 nm to 1000 nm in organic matter and clay minerals. The right peak is very high, corresponding to a pore size of 40-200 µm. The highly developed lamellation in shale formations generates micro fractures, which may happen during sample preparation or mercury injection tests. Therefore, there is a high probability that the right peak corresponds to artificial fractures. Considering the fact that large pores (>5 µm) correspond to mercury injection pressure less than 0.14 MPa, these artificial fractures are more likely to be induced during sample preparation. Figure 12 shows that macropores account for the largest percentage (73.17%), Relative content, % Pore diameter, nm X3-17 is advantageous for the storage of adsorbed gas. Meanwhile, we see that pores larger than 10 nm are not very developed. In a PSD frequency histogram ( Figure 10), mesopores account for the largest percentage of all pores, followed by micropores and macropores. Micropores and mesopores accounted for more than 90% of the total pore volume.

Mercury Intrusion Results
Two peaks can be found on the MICP measurement results of PSD, as shown in Figure 11. The left peak is relatively small and smooth, corresponding to macropores of 10 nm to 1000 nm in organic matter and clay minerals. The right peak is very high, corresponding to a pore size of 40-200 µm. The highly developed lamellation in shale formations generates micro fractures, which may happen during sample preparation or mercury injection tests. Therefore, there is a high probability that the right peak corresponds to artificial fractures. Considering the fact that large pores (>5 µm) correspond to mercury injection pressure less than 0.14 MPa, these artificial fractures are more likely to be induced during sample preparation. Figure 12 shows that macropores account for the largest percentage (73.17%), Relative content, % Pore diameter, nm X3-17

Mercury Intrusion Results
Two peaks can be found on the MICP measurement results of PSD, as shown in Figure 11. The left peak is relatively small and smooth, corresponding to macropores of 10 nm to 1000 nm in organic matter and clay minerals. The right peak is very high, corresponding to a pore size of 40-200 µm. The highly developed lamellation in shale formations generates micro fractures, which may happen during sample preparation or mercury injection tests. Therefore, there is a high probability that the right peak corresponds to artificial fractures. Considering the fact that large pores (>5 µm) correspond to mercury injection pressure less than 0.14 MPa, these artificial fractures are more likely to be induced during sample preparation. Figure 12 shows that macropores account for the largest percentage (73.17%), followed by mesopores (26.83%). Due to the limitations of the instrument, no micropores are detected by MICP. Whether the tested macropores are primitive or induced needs to be determined by combining with other techniques, such as NMR.

Nuclear Magnetic Resonance Results
Two or three peaks can be found on PSD, measured by NMR ( Figure 13). The left peak, corresponding to pores smaller than 10 nm, has the largest percentage, while the other two peaks correspond to larger pore sizes of 800 nm and 7000 nm, respectively. The NMR results indicate that small pores are very developed in shale formations, while large pores account for a small but nonnegligible percentage ( Figure 14). Generally, the left two peaks correspond to micro-, meso-, and macro-pores in matrix, while the third peak corresponds to micro fractures. The PSD histogram measured by the NMR is similar to that of N2 adsorption in terms of the small pore size range, while large pores could not be detected by N2 adsorption, but by NMR.

Nuclear Magnetic Resonance Results
Two or three peaks can be found on PSD, measured by NMR ( Figure 13). The left peak, corresponding to pores smaller than 10 nm, has the largest percentage, while the other two peaks correspond to larger pore sizes of 800 nm and 7000 nm, respectively. The NMR results indicate that small pores are very developed in shale formations, while large pores account for a small but non-negligible percentage ( Figure 14). Generally, the left two peaks correspond to micro-, meso-, and macro-pores in matrix, while the third peak corresponds to micro fractures. The PSD histogram measured by the NMR is similar to that of N 2 adsorption in terms of the small pore size range, while large pores could not be detected by N 2 adsorption, but by NMR.

Comprehensive Analysis of Pore Size Distribution
N2 adsorption, MICP, and NMR can all measure PSD and reflect the heterogeneity of shale samples, with different test ranges. The lowest limit of N2 adsorption is 0.35 nm, while the MICP test range is 3.6 nm-950 µm, and the NMR test range is 1 nm-5 mm. Comparing the results from the three methods, we find that the N2 adsorption results mainly reflect the micropore and mesopore size distribution, while the NMR results reflect all pore size range and MICP results mainly test the development of macropores and micro fractures. Although the N2 adsorption and NMR results display similar PSD trends of small pore size ranges, the peak of N2 adsorption results (3-4 nm) is slightly smaller than that of NMR (4-5 nm). This is because samples are saturated by plant oil in the NMR test, and they struggle to enter micropores due to their large diameter compared to nitrogen molecules. Therefore, the micropore size determined by NMR is larger than that measured by N2 adsorption. The nitrogen molecule is much smaller, so it enters micropores more easily than oil molecules. Therefore, the N2 adsorption results are closer to the real data compared to NMR results.

Comprehensive Analysis of Pore Size Distribution
N 2 adsorption, MICP, and NMR can all measure PSD and reflect the heterogeneity of shale samples, with different test ranges. The lowest limit of N 2 adsorption is 0.35 nm, while the MICP test range is 3.6 nm-950 µm, and the NMR test range is 1 nm-5 mm. Comparing the results from the three methods, we find that the N 2 adsorption results mainly reflect the micropore and mesopore size distribution, while the NMR results reflect all pore size range and MICP results mainly test the development of macropores and micro fractures. Although the N 2 adsorption and NMR results display similar PSD trends of small pore size ranges, the peak of N 2 adsorption results (3-4 nm) is slightly smaller than that of NMR (4-5 nm). This is because samples are saturated by plant oil in the NMR test, and they struggle to enter micropores due to their large diameter compared to nitrogen molecules. Therefore, the micropore size determined by NMR is larger than that measured by N 2 adsorption. The nitrogen molecule is much smaller, so it enters micropores more easily than oil molecules. Therefore, the N 2 adsorption results are closer to the real data compared to NMR results. N 2 adsorption can measure micropores and mesopores accurately, while MICP mainly tests macropores. Consequently, combining the two methods can better characterize PSD in shale formations.  N2 adsorption can measure micropores and mesopores accurately, while MICP mainly tests macropores. Consequently, combining the two methods can better characterize PSD in shale formations. Figures 15 and 16 show the PSD of the Longmaxi shale samples, tested comprehensively by N2 adsorption and MICP, where pores smaller than 50 nm are measured by N2 adsorption and pores larger than 50 nm are measured by MICP. Pore volume in the Longmaxi Shale Formation is mainly mesopores and macropores (including artificial fractures).

Gas Composition and Origin
Milkov et al. [37] studied gas composition and origins based on around 2600 shale gas samples from 76 geological formations in 38 sedimentary basins located in eleven countries. It is found that methane is the predominated hydrocarbon component, with more than 80% in volume concentration, followed by ethane with around 6% in volume concentration, and propane with around 2% in volume

Gas Composition and Origin
Milkov et al. [37] studied gas composition and origins based on around 2600 shale gas samples from 76 geological formations in 38 sedimentary basins located in eleven countries. It is found that methane is the predominated hydrocarbon component, with more than 80% in volume concentration, followed by ethane with around 6% in volume concentration, and propane with around 2% in volume concentration. Nitrogen and carbon dioxide are two main non-hydrocarbon components in shale gas samples, with average volume concentration around 6% and 2%, respectively. For most shale plays in the USA, China and Argentina, it is found that the most productive and commercially successful shale Energies 2020, 13, 5427 16 of 50 plays have pure thermogenic origin. This is very different from the study of Curtis [19], where it is found that shale gas has predominantly microbial origin.

Shale Gas Occurrence Types
The types of natural gas in shale formations are determined by diverse formation physics and pore characteristics. In accordance with the classification of pore structures in shale formations (Section 2.3), free gas, adsorbed gas, and dissolved gas are three possible gas occurrence states underground [38,39]. Generally, free gas is stored not only in fractures, but also in pore systems, including organic pores and inorganic pores. Adsorbed gas is mostly stored on the surface of organic matter in equilibrium state with free gas. Dissolved gas is usually stored in liquid hydrocarbons, formation water, but most importantly in solid kerogen. Organic kerogen serves as the source rock and generates shale gas continuously [40].
The percentage of different gas types varies from one reservoir to another, since it is significantly influenced by pressure, temperature, organic matter types, organic matter content, organic matter maturity, the development of micro fractures, and liquid hydrocarbon content. The different gas types and gas flow mechanisms in organic shale nanopores can be seen in Figure 17. concentration. Nitrogen and carbon dioxide are two main non-hydrocarbon components in shale gas samples, with average volume concentration around 6% and 2%, respectively. For most shale plays in the USA, China and Argentina, it is found that the most productive and commercially successful shale plays have pure thermogenic origin. This is very different from the study of Curtis [19], where it is found that shale gas has predominantly microbial origin.

Shale Gas Occurrence Types
The types of natural gas in shale formations are determined by diverse formation physics and pore characteristics. In accordance with the classification of pore structures in shale formations (Section 2.3), free gas, adsorbed gas, and dissolved gas are three possible gas occurrence states underground [38,39]. Generally, free gas is stored not only in fractures, but also in pore systems, including organic pores and inorganic pores. Adsorbed gas is mostly stored on the surface of organic matter in equilibrium state with free gas. Dissolved gas is usually stored in liquid hydrocarbons, formation water, but most importantly in solid kerogen. Organic kerogen serves as the source rock and generates shale gas continuously [40].
The percentage of different gas types varies from one reservoir to another, since it is significantly influenced by pressure, temperature, organic matter types, organic matter content, organic matter maturity, the development of micro fractures, and liquid hydrocarbon content. The different gas types and gas flow mechanisms in organic shale nanopores can be seen in Figure 17.

Free Gas Characterization
Free gas is stored in organic or inorganic pores, micro fractures, and hydraulic fractures. The content of free gas is determined by adsorbed gas and dissolved gas. Only when the total gas amount is larger than the sum of adsorbed and dissolved gas amount does a free gas state exist. Under highpressure or high-temperature reservoir conditions, gas behaviors do not satisfy the ideal gas equation of state (EOS), so the real gas EOS needs to be adopted to describe its behaviors [42,43].

Semi-Empirical Formula
EOS is one of the most important models to calculate the thermodynamic properties of real gas. More than 150 types of EOS have been proposed to describe the pressure-volume-temperature relationships of real gas, but none is able to include all properties of any gas in the engineering application range [44]. Typical and widely applied EOS, as well as their pros and cons, are displayed in Table 4.

Free Gas Characterization
Free gas is stored in organic or inorganic pores, micro fractures, and hydraulic fractures. The content of free gas is determined by adsorbed gas and dissolved gas. Only when the total gas amount is larger than the sum of adsorbed and dissolved gas amount does a free gas state exist. Under high-pressure or high-temperature reservoir conditions, gas behaviors do not satisfy the ideal gas equation of state (EOS), so the real gas EOS needs to be adopted to describe its behaviors [42,43].

Semi-Empirical Formula
EOS is one of the most important models to calculate the thermodynamic properties of real gas. More than 150 types of EOS have been proposed to describe the pressure-volume-temperature relationships of real gas, but none is able to include all properties of any gas in the engineering application range [44]. Typical and widely applied EOS, as well as their pros and cons, are displayed in Table 4.

Empirical Formula
According to the corresponding state principle, gas deviation factor Z could be introduced to describe real gas behaviors: where p is the pressure, Pa; V is the gas volume, m 3 ; Z is the gas deviation factor, dimensionless; n is the molar mass, kg/mol; R is the universal gas constant, J/(mol K); and T is temperature, K.
In Equation (8), gas deviation factor Z could be obtained by experiments, referring to Z-plate [44], or calculating from an empirical formula [27]. The Z-plate and empirical formula were mainly obtained by assuming Z c as a constant in the range of 0.23-0.29. That is to say, gas deviation factor Z is a function of the reduced pressure p r and reduced temperature T r . Therefore, we also name the Z-plate as the two-parameter generalized compressibility chart. The critical gas deviation factor Z c of most materials varies in the range of 0.23-0.29. Therefore, a more precise expression of Z is expected to be obtained by regarding Z as a function of p r , T r , and another parameter (Z c or acentric factor ω)-a three-parameter relationship [44].

Adsorbed Gas
Adsorbed gas is mainly stored on the surface of matrix particles, kerogen, and clay minerals, and can account for 20-85% of total gas reserves [19,[45][46][47]. Gas adsorption on shale matrix particles belongs to physical adsorption [48]. Although adsorbed gas contributes to the total gas production, its exact contribution is not clear. Compared to the contribution of adsorbed gas in total gas production, the percentage of adsorbed gas in original gas in place (OGIP) is much clearer. Recent studies have found that adsorbed gas accounts for 50-80% of OGIP when the pressure is lower than 13.79 MPa, while it accounts for 30-50% when the pressure is higher than 13.79 MPa [49,50].
Organic nanopores of 5-750 nm, are significantly developed in shale gas reservoirs when the maturity of organic matter is larger than 0.6% [28,30]. Due to the small pore radius, organic matter has a large surface area. For example, the specific surface area can be up to 300 m 2 /g in nanoporous kerogen [51]. The enormous surface area provides favorable places for gas adsorption. The gas adsorption capacity is mainly affected by pore structures, mineral compositions, metamorphism degree, gas components, pressure, temperature, water vapor content, etc. Since adsorbed gas is mainly stored in organic matter, the TOC content significantly affects the adsorbed gas content in shale gas reservoirs. As we can see from Figure 18, there is a positive relationship between the adsorbed gas content and TOC content in different shale gas reservoirs. This is because large TOC content means more organic matter in the shale matrix, which can provide sufficient space for gas storage due to its large surface area.
Taking the Barnett Shale as an example, as shown in Figure 18, we can distinguish different gas types from the relationship between adsorbed gas and the total gas amount in shale matrix. Adsorbed gas and free gas in the organic matrix increase as the TOC content increases, while the amount of free gas in the inorganic matrix is not affected by the TOC content. In shale gas plays, adsorbed gas is a non-ignorable component in shale gas reserve calculation [27,39,[53][54][55], and gas ad-/desorption is important in the study of gas flow behaviors [56][57][58][59]. If we analyze the case further, we could ask how much adsorbed gas could be produced during shale gas production, and how significantly gas ad-/desorption affects gas transient flow behaviors in shale gas reservoirs. Taking the Barnett Shale as an example, as shown in Figure 18, we can distinguish different gas types from the relationship between adsorbed gas and the total gas amount in shale matrix. Adsorbed gas and free gas in the organic matrix increase as the TOC content increases, while the amount of free gas in the inorganic matrix is not affected by the TOC content. In shale gas plays, adsorbed gas is a non-ignorable component in shale gas reserve calculation [27,39,[53][54][55], and gas ad-/desorption is important in the study of gas flow behaviors [56][57][58][59]. If we analyze the case further, we could ask how much adsorbed gas could be produced during shale gas production, and how significantly gas ad-/desorption affects gas transient flow behaviors in shale gas reservoirs.
The gas transporting process of CH4 and He in organic shale samples was compared by an experimental study at 3.4 MPa and 308 K [60], where CH4 serves as the adsorptive gas and He is nonadoptive. The production dynamics of CH4 and He can be seen in Figure 19a. Assuming the free gas amount of CH4 and He is equal in shale samples, the adsorbed volume of CH4 can be obtained by the difference between the total CH4 production volume and the total He production volume, which are 2.60 cm 3 /g and 1.33 cm 3 /g respectively. Therefore, the produced volume of free gas for unit mass shale particles under standard conditions is 1.33 cm 3 /g, while the adsorbed gas amount is 1.27 cm 3 /g. Similarly, simulation results from dynamic adsorption diffusion model show that the production of free gas dominates at an early production period (before point A) and drops very fast, while adsorbed gas dominates the later production after point A for a relatively long time. Experimental study and model simulation signified that both free gas and adsorbed gas played an important role in gas production.
The above research is conducted at low pressure (3.4 MPa) and temperature (308 K) compared to practical shale gas reservoir conditions. Gas desorption pressure in shale is usually below 12 MPa, which is close to the abandonment pressure of shale gas reservoirs. Meanwhile, the formation pressure is mainly depleted in a small area near the wellbore, i.e., the average pressure of shale formation is much higher than the abandon pressure. Consequently, gas desorption may only occur in a small area near the wellbore or hydraulic fractures, meaning a limited amount of adsorbed gas is produced during the life cycle of the shale gas reservoir. The significance of adsorbed gas, as well as the corresponding seepage mechanisms, need further investigation. The gas transporting process of CH 4 and He in organic shale samples was compared by an experimental study at 3.4 MPa and 308 K [60], where CH 4 serves as the adsorptive gas and He is non-adoptive. The production dynamics of CH 4 and He can be seen in Figure 19a. Assuming the free gas amount of CH 4 and He is equal in shale samples, the adsorbed volume of CH 4 can be obtained by the difference between the total CH 4 production volume and the total He production volume, which are 2.60 cm 3 /g and 1.33 cm 3 /g respectively. Therefore, the produced volume of free gas for unit mass shale particles under standard conditions is 1.33 cm 3 /g, while the adsorbed gas amount is 1.27 cm 3 /g. Similarly, simulation results from dynamic adsorption diffusion model show that the production of free gas dominates at an early production period (before point A) and drops very fast, while adsorbed gas dominates the later production after point A for a relatively long time. Experimental study and model simulation signified that both free gas and adsorbed gas played an important role in gas production.
The above research is conducted at low pressure (3.4 MPa) and temperature (308 K) compared to practical shale gas reservoir conditions. Gas desorption pressure in shale is usually below 12 MPa, which is close to the abandonment pressure of shale gas reservoirs. Meanwhile, the formation pressure is mainly depleted in a small area near the wellbore, i.e., the average pressure of shale formation is much higher than the abandon pressure. Consequently, gas desorption may only occur in a small area near the wellbore or hydraulic fractures, meaning a limited amount of adsorbed gas is produced during the life cycle of the shale gas reservoir. The significance of adsorbed gas, as well as the corresponding seepage mechanisms, need further investigation.
Assuming adsorbed gas volume is a function of pressure, Tang et al. [61] obtained the absolute adsorbed gas amount from excess adsorption and studied the adsorbed gas proportion to total gas at different shale depths ( Figure 20). The conventional absolute adsorption refers to results obtained by fitting low and intermediate pressure sorption data using the Langmuir model of Equation (10), while the new absolute adsorption refers to a dual-site Langmuir model considering the adsorbed layer variation and excess adsorption. The conventional model severely underestimates the absolute adsorption amount when the pressure is higher than 6 MPa, as shown in Figure 20a. The percentage of adsorbed gas to total gas in place (GIP) is a function of shale formation depth, where it increases fast in shallow areas and slows down after 2000 m, as shown in Figure 20b. The adsorbed gas accounts for approximately 40-80% at different depths of formation. Meanwhile, the excess adsorption amount needs to be corrected to the absolute adsorption amount when considering the adsorbed gas percentage Assuming adsorbed gas volume is a function of pressure, Tang et al. [61] obtained the absolute adsorbed gas amount from excess adsorption and studied the adsorbed gas proportion to total gas at different shale depths ( Figure 20). The conventional absolute adsorption refers to results obtained by fitting low and intermediate pressure sorption data using the Langmuir model of Equation (10), while the new absolute adsorption refers to a dual-site Langmuir model considering the adsorbed layer variation and excess adsorption. The conventional model severely underestimates the absolute adsorption amount when the pressure is higher than 6 MPa, as shown in Figure 20a. The percentage of adsorbed gas to total gas in place (GIP) is a function of shale formation depth, where it increases fast in shallow areas and slows down after 2000 m, as shown in Figure 20b. The adsorbed gas accounts for approximately 40-80% at different depths of formation. Meanwhile, the excess adsorption amount needs to be corrected to the absolute adsorption amount when considering the adsorbed gas percentage in GIP. Otherwise, it will massively underestimate the adsorbed gas amount and overestimate the free gas amount.

Dissolved Gas
After the equilibrium between adsorption and desorption is found, shale gas could dissolve into the liquid hydrocarbon or formation water during the hydrocarbon accumulation process. Meanwhile, organic kerogen continuously generates shale gas and contains a certain amount of gas molecules [62]. The gas in liquid hydrocarbon, formation water, and kerogen is called the dissolved gas, which has been overlooked, but may play a significant role in shale gas reservoir development [40,59]. Gap-filling and hydration are the two main storage mechanism of dissolved gas, and can be described theoretically by Henry's law [63]: where C b is the mole concentration of dissolved gas, mol/m 3 and K c is the Henry constant, m 3 Pa/mol. Since it is hard to differentiate dissolved gas from adsorbed gas, both gas types are usually attributed to one type, namely adsorbed gas. Moreover, adsorbed gas and dissolved gas can be transformed to the other under proper circumstances. Therefore, it can be roughly seen as one type in some cases [48,64].

Dissolved Gas
After the equilibrium between adsorption and desorption is found, shale gas could dissolve into the liquid hydrocarbon or formation water during the hydrocarbon accumulation process. Meanwhile, organic kerogen continuously generates shale gas and contains a certain amount of gas molecules [62]. The gas in liquid hydrocarbon, formation water, and kerogen is called the dissolved gas, which has been overlooked, but may play a significant role in shale gas reservoir development [40,59]. Gap-filling and hydration are the two main storage mechanism of dissolved gas, and can be described theoretically by Henry's law [63]: where Cb is the mole concentration of dissolved gas, mol/m 3 and Kc is the Henry constant, m 3 Pa/mol.

Gas Adsorption and Desorption
Shale gas can be stored on pore surfaces of organic matter and clays by gas adsorption. Organic matter in the shale matrix is a key parameter that influences gas adsorption characteristics in shale gas reservoirs. On the one hand, a large amount of nanopores are developed in organic shales, which provide enormous surface area for the gas to be adsorbed on. On the other hand, the adsorption potential is significant in organic nanopores compared to in inorganic nanopores or large organic pores. The adsorption-desorption law in organic shale nanopores is a key scientific problem in the practice of shale gas development, affecting the accuracy of evaluating shale adsorption capacity, studying the seepage flow behaviors, and developing transient seepage mathematical models [65].

Different Sorption Types and Models
Methane is the main component of shale gas underground, with a critical pressure of 4.59 MPa and a critical temperature of 190.53 K. Therefore, shale gas is in a supercritical state under in situ formation conditions (3000-6000 m, with high pressure up to 60 MPa) [66]. The study of supercritical gas sorption is essential for an accurate understanding of adsorption and desorption mechanisms in shale gas reservoirs. Gas sorption mechanisms are quite confusing, and no unified conclusions have been reached. Monolayer adsorption, multilayer adsorption, and micropore filling are three common assumptions in shale gas sorption research. Based on these assumptions, the Langmuir model, BET model, the Dubinin-Radushkevich (D-R) model, and the Dubinin-Astakhov (D-A) model have been established to fit the sorption data, and have obtained good results. However, good fitting results do not guarantee the validity of the assumption in the adsorption models. For example, even if the Langmuir model fits the experimental data very well, we cannot say that gas adsorption belongs to monolayer adsorption.

Monolayer Adsorption Type
Monolayer adsorption means that gas molecules adsorb on the pore surface in one layer, so the thickness of adsorbed gas equals the molecular diameter. Due to the huge surface area in a shale matrix, a considerable amount of adsorbed gas exists in shale gas reservoirs. Assuming 80% gas saturation in shale samples, the adsorbed gas ratio can be 22.65% of the total gas amount according to the Ono-Kondo lattice model established by Zhou et al. [67].
Due to its simplicity and accurate fitting to the experimental data, it is widely used to describe monolayer gas adsorption, which can be written as: where G a is the absolute adsorbed gas amount; G m is the maximum adsorbed capacity, p is the pressure, and b is the Langmuir sorption constant, which can be obtained by fitting the experimental or field test data.
(2) Freundlich model [69]: With pressure decreasing, the Langmuir equation of Equation (10) approaches the Henry's law of Equation (9). Therefore, the Henry's law can describe low-pressure sorption behaviors, since any sorption isotherm satisfies the linear relationship between adsorption amount and pressure at low pressure. To broaden the application range of the Henry's law into high-pressure areas, an exponential empirical formulation, namely the Freundlich model, was used: where k f is related to the adsorption interaction and adsorption amount; n f is a constant usually between 2 and 3, reflecting the intensity of adsorption. The values of both k f and n f depend on the type of adsorbent and adsorbate as well as the temperature.
With temperature increasing, constant n f approaches to unity and the Freundlich model of Equation (11) becomes the Henry model of Equation (9). The Freundlich model can properly describe monolayer adsorption, especially for low-concentration gases and in the meso pressure range. However, there is no explicit physical meaning of the constants k f and n f , and it cannot explain the mechanisms of adsorption.
where l reflects the heterogeneity of adsorbents, l ≤ 1. The smaller the value of l, the stronger the heterogeneity of the adsorbent. If an adsorbent possesses an ideal surface, l tends to be 1 and the Langmuir-Freundlich equation is equivalent to the Langmuir equation. (4) Toth model [69]: To improve the fitting capacity of the Langmuir model, the Toth model was proposed: where t is a constant related to the adsorbent properties.
Note that the Toth model of Equation (13) solves two problems: (1) the Freundlich model of Equation (11) and the DR model of Equation (16) do not satisfy Henry's law at low pressure; and (2) no maximum adsorption amount appears in the Freundlich model of Equation (11) with increasing pressure. Bae et al. [69] found that the Toth equation fitted the experimental data better than the extended three-parameter and Langmuir equation, and yielded realistic values of pore volumes of coal samples and adsorbed gas density.

Multilayer Sorption Type
(1) Two-parameter BET model [32]: Multilayer adsorption can be modeled by the BET sorption theory, which assumes that gas molecules can adsorb on a solid surface by infinite layers and no interaction exists between contiguous layers. In other words, any monolayer obeys the Langmuir adsorption theory in the BET model, which can be expressed as follows: where G m is the maximum monolayer adsorption amount, p s is the saturated vapor pressure, and C b is the dimensionless constant controlling the time of multilayer adsorption.
In applying Equation (14), its equivalent expression needs to be adopted, which is Equation (1). A plot of (p/p s )/[G a (1 − p/p s )] versus p/p s is employed to figure out whether the adsorption follows the BET theory. If the plot satisfies the linear relationship in the range of 0.005 < p/p s < 0.35, we can use the scope and intercept of the straight line to obtain the values of G m and C.
(2) Three-parameter BET model [44,69]: The above two-parameter BET model of Equation (14) assumes infinite layers of adsorption. If the adsorption layers are finite, the three-parameter BET model can be employed: If n b = 1, Equation (15) simplifies into the Langmuir model of Equation (10); if n→∞, Equation (15) transforms into the two-parameter BET model of Equation (14). Note that the three-parameter BET model of Equation (15) is applicable to describe adsorption behaviors for p/p s in the range of 0.35-0.60.
Note that, although in theory the multilayer adsorption assumption may produce a wider application scope with the BET model than the Langmuir model, it may not be suitable to adopt the Energies 2020, 13, 5427 23 of 50 BET theory in a shale gas adsorption study, since sorption behaviors in shale gas reservoirs belong to supercritical adsorption and the saturated vapor pressure p s of shale gas (mainly methane) does not exist in practice. Besides, as reported, the BET model may have a poorer performance than the Langmuir model at fitting absolute sorption data, as reported by [70].

Dubinin-Radushkevich and Dubinin-Astakhov Models
Gas molecule behavior in nanopores is significantly different from that in mesopores or macropores, since there is a superposition of adsorption potential from both pore sides. Consequently, the adsorption force of micropore walls on gas molecules is much greater than in mesopores or macropores, leading to large adsorption. Dubinin named the gas adsorption in these small-scale pores micropore volume filling [71]. Compared to the Langmuir adsorption theory, micropore volume filling is more helpful for understanding the gas adsorption mechanism and gas true storage forms, and for evaluating gas adsorption properties. Meanwhile, it has been reported that the D-A model provides a better fit to sorption data of coal than the Langmuir model [72].
The condensed adsorbate looks like microemulsion droplets when adsorption occurs in micropores, which is greatly affected by interfaces. Based on the Polanyi adsorption potential theory [73], the D-R model and D-A model are commonly used in shale gas adsorption studies, and can be expressed as follows: where W is the pore volume filled with gas molecules at relative pressure p/p s ; W 0 is the total volume of micropores; D is a parameter related to the adsorbate-adsorbent system; and m is a parameter ranging from 2 to 6, reflecting the heterogeneity of potential energy on adsorbent surfaces. If m constantly equals 2, Equation (16) is the D-R model. If m is a random parameter between 2 and 6, then Equation (16) is the D-A model. In applying Equation (16), W is equivalent to the absolute adsorbed gas amount G a , and W 0 is equivalent to the maximum adsorbed gas amount G m , Equation (16) can be re-expressed as follows: Compared to the D-R model, the D-A model performs better fitting with experimental data, according to the study of Wang et al. [70]. This is because the chosen range of structure heterogeneity parameter m in the D-A model is broader than that in the D-R model, which is related to pore size distribution in shale formations. The chosen parameter m in the D-A model brings the micropore structure information into adsorption prediction and modeling, while it has a constant value of 2 in the D-R model, without considering the structural heterogeneity in shale samples.

Calculation of Virtual Saturated Vapor Pressure
Note that the above D-R and D-A models were not initially proposed for supercritical adsorption, but for subcritical sorption. Therefore, the concept of saturated vapor pressure in the D-R and D-A model was replaced by virtual saturated vapor pressure or supercritical adsorption limited pressure [72]. Generally, the virtual saturated vapor pressure can be calculated by the following empirical formulations or approaches: (1) The first is the Dubinin method [71,74,75]: where p c is the critical pressure and T c is the critical temperature. (2) The second is the Reid method [72]: where T b is the boiling point of gas at atmospheric pressure. (3) The third is the Antoine method [74]: Equation (20) where parameters c A and d A are determined by the gas critical point (T c , p c ) and boiling point (T b , 101,325 Pa). For methane, c A = −1032.693 and d A = 6.945. (5) The Amankwah method [76] is an improved calculation method for the Dubinin method, which involves a parameter k A to account for interactions in an adsorbate-adsorbent system: where k A is a parameter accounting for interactions in the adsorbate-adsorbent systems. (6) The linearization of isotherm adsorption data is another processing method [77] to extend the D-R and D-A models into supercritical area. By transforming isotherms from G ex versus p space to ln[lnG ex ] − 1 versus lnp space, where G ex is the excess adsorption amount, a bunch of fitting straight lines could be obtained; they converge to a single point B, as shown in Figure 21. This merge point B is defined as the limiting state of the adsorbate, corresponding to the extreme condition of the adsorption potential field, where no more adsorptive molecules can enter the adsorbent micropores [77]. Therefore, the limiting pressure and limiting adsorption amount of the merge point correspond to the saturated vapor pressure and the saturated adsorption amount in the D-R and D-A models, respectively.
In order to compare these six different methods, the experimental data in Zhang's study for organic-rich Woodford shale [78] are collected for analysis at temperatures of 308.53 K, 323.53 K, and 338.53 K, respectively, as shown in Figure 22a. The linearization processing method of adsorption data in Zhou's study is adopted, which transforms the sorption data into three straight lines and defines a limiting state at the intersection point A in Figure 22b. The other five methods are also employed to calculate the virtual saturated vapor pressure for the same shale sample and sorption data. The calculated results are shown in Figure 23, which exhibits great differences to the results from different methods. The results from Antoine and Astakhov have obviously larger values than the others, while the linearization processing results, which are independent of temperature, have much smaller values than the others. Dubinin and Reid's results are basically the same, with slightly lower values than in the Amankwah results. Dubinin's method only considers the properties of the individual adsorbates, while the Amankwah method also takes adsorbent properties into consideration. Moreover, the value of parameter k is obtained by nonlinear fitting for isotherm sorption data, which is more practical than the constant value in the Dubinin method. As a result, the Amankwah method is recommended to calculate virtual saturated vapor pressure.
In order to compare these six different methods, the experimental data in Zhang's study for organic-rich Woodford shale [78] are collected for analysis at temperatures of 308.53 K, 323.53 K, and 338.53 K, respectively, as shown in Figure 22a. The linearization processing method of adsorption data in Zhou's study is adopted, which transforms the sorption data into three straight lines and defines a limiting state at the intersection point A in Figure 22b. The other five methods are also employed to calculate the virtual saturated vapor pressure for the same shale sample and sorption data. The calculated results are shown in Figure 23, which exhibits great differences to the results from different methods. The results from Antoine and Astakhov have obviously larger values than the others, while the linearization processing results, which are independent of temperature, have much smaller values than the others. Dubinin and Reid's results are basically the same, with slightly lower values than in the Amankwah results. Dubinin's method only considers the properties of the individual adsorbates, while the Amankwah method also takes adsorbent properties into consideration. Moreover, the value of parameter k is obtained by nonlinear fitting for isotherm sorption data, which is more practical than the constant value in the Dubinin method. As a result, the Amankwah method is recommended to calculate virtual saturated vapor pressure.
Note that the linearization processing method was proposed to tackle sorption problems in the ranges of 77-298 K and 0-7 MPa [77]. The pressure and temperatures in shale gas reservoirs are generally beyond this range, so this method is not recommended in shale gas sorption studies.  T=77K  T=93K  T=113K   T=133K  T=153K  T=173K   T=193K  T=213K  T=233K   T=253K  T=273K  . Linearized adsorption isotherms of hydrogen [77], where p is in KPa and G ex is in mmol/g. Figure 21. Linearized adsorption isotherms of hydrogen [77], where p is in KPa and Gex is in mmol/g.
In order to compare these six different methods, the experimental data in Zhang's study for organic-rich Woodford shale [78] are collected for analysis at temperatures of 308.53 K, 323.53 K, and 338.53 K, respectively, as shown in Figure 22a. The linearization processing method of adsorption data in Zhou's study is adopted, which transforms the sorption data into three straight lines and defines a limiting state at the intersection point A in Figure 22b. The other five methods are also employed to calculate the virtual saturated vapor pressure for the same shale sample and sorption data. The calculated results are shown in Figure 23, which exhibits great differences to the results from different methods. The results from Antoine and Astakhov have obviously larger values than the others, while the linearization processing results, which are independent of temperature, have much smaller values than the others. Dubinin and Reid's results are basically the same, with slightly lower values than in the Amankwah results. Dubinin's method only considers the properties of the individual adsorbates, while the Amankwah method also takes adsorbent properties into consideration. Moreover, the value of parameter k is obtained by nonlinear fitting for isotherm sorption data, which is more practical than the constant value in the Dubinin method. As a result, the Amankwah method is recommended to calculate virtual saturated vapor pressure.
Note that the linearization processing method was proposed to tackle sorption problems in the ranges of 77-298 K and 0-7 MPa [77]. The pressure and temperatures in shale gas reservoirs are generally beyond this range, so this method is not recommended in shale gas sorption studies.

Gibbs Excess Adsorption
For high-pressure and -temperature sorption, Gibbs excess adsorption is adopted to describe its unique behaviors [79]. Generally, adsorbed gas and bulk gas both exist for an adsorbate-adsorbent system, where adsorbed gas is distributed on the pore surface as a layer and bulk gas is far from the surface. Bulk phase gas is also distributed in the adsorption layer, which is irrelevant to gas-solid molecular interactions and can be ignored at low pressure. However, it needs to be considered in shale gas sorption research, since the in situ pressure is high in shale gas reservoirs (>30 MPa) [65]. Therefore, the excess adsorption amount corresponds to the part that is larger than the bulk phase density in the Note that the linearization processing method was proposed to tackle sorption problems in the ranges of 77-298 K and 0-7 MPa [77]. The pressure and temperatures in shale gas reservoirs are generally beyond this range, so this method is not recommended in shale gas sorption studies.

Gibbs Excess Adsorption
For high-pressure and -temperature sorption, Gibbs excess adsorption is adopted to describe its unique behaviors [79]. Generally, adsorbed gas and bulk gas both exist for an adsorbate-adsorbent system, where adsorbed gas is distributed on the pore surface as a layer and bulk gas is far from the surface. Bulk phase gas is also distributed in the adsorption layer, which is irrelevant to gas-solid molecular interactions and can be ignored at low pressure. However, it needs to be considered in shale gas sorption research, since the in situ pressure is high in shale gas reservoirs (>30 MPa) [65]. Therefore, the excess adsorption amount corresponds to the part that is larger than the bulk phase density in the adsorption layer. The difference between the absolute adsorption and the excess adsorption can be seen in Figure 25, where the absolute adsorption (Figure 25e) consists of excess adsorption (Figure 25c) and bulk phase gas in the adsorbed layer (Figure 25d). The relationship can be expressed as follows: where G ex is the Gibbs excess adsorption amount, G a is the absolute adsorption amount, and v ad is the adsorbed gas volume, as can be seen in Figure 25. The measured adsorption amount in high-pressure sorption experiments is the Gibbs excess adsorption amount G ex in Equation (23). The relationship has another explanation: namely, the adsorbed gas is under the effect of bulk phase gas buoyancy. Thus, the measured adsorbed gas weight equals the difference between absolute adsorbed gas weight and the buoyancy it received in bulk-phase gas.
For adsorbed gas, the following relationship exists: where ρ ad is the adsorbed gas density. Then, incorporating Equations (23) and (24), we can obtain: . (25) Associating Equation (25) with the calculation methods for bulk and adsorbed gas density in Sections 3.4.1 and 3.4.2, the simulation results of absolute adsorption amount G a can be transformed into measured excess adsorption amounts G ex . In low-pressure sorption studies, the bulk gas density ρ g is much lower than the adsorbed gas density ρ ad , and the excess adsorption amount G ex is approximately the same as the absolute adsorption amount G a . However, gas density ρ g becomes comparable to adsorbed gas density ρ ad with pressure increasing, as to be introduced in Section 3.4.2. Thus, the difference between the absolute adsorption amount G a and the excess adsorption amount G ex cannot be ignored and there is a maximum on the plot of measured adsorption amount versus pressure or bulk phase gas density, as can be seen from Figure 24a. The location of the maximum depends on the interaction between adsorbate and adsorbent, as well as the thermodynamic state of the adsorptive [80].

Supercritical Adsorption Models
The introduction of the virtual saturated vapor pressure in Section 3.1.3 extends the micropore filling models from subcritical area into supercritical range. However, Sakurovs et al. [81] noted that this method cannot easily accommodate adsorption at conditions where both the pressure and  Figure 24. Fitting of models to CO 2 adsorption data (a) and residuals of curves from different models to CH 4 adsorption data, and the base line means the measured equals to the predicted (b) [81].
approximately the same as the absolute adsorption amount Ga. However, gas density ρg becomes comparable to adsorbed gas density ρad with pressure increasing,as to be introduced in Section 3.4.2. Thus, the difference between the absolute adsorption amount Ga and the excess adsorption amount Gex cannot be ignored and there is a maximum on the plot of measured adsorption amount versus pressure or bulk phase gas density, as can be seen from Figure 25a. The location of the maximum depends on the interaction between adsorbate and adsorbent, as well as the thermodynamic state of the adsorptive [80].

Supercritical Adsorption Models
The introduction of the virtual saturated vapor pressure in Section 3.1.3 extends the micropore filling models from subcritical area into supercritical range. However, Sakurovs et al. [81] noted that this method cannot easily accommodate adsorption at conditions where both the pressure and temperature are above the critical values. Since the adsorbed gas density is greater than the free gas density in supercritical conditions, another method, which replaced the saturated vapor pressure p s by adsorbed phase gas density ρ ad , and gas pressure p by gas density ρ g , was proposed to extend the volume-filling models to a wider pressure and temperature application range [72]. Based on this idea and associating Equation (17) with Equation (25), we obtain: Similarly, other models in Section 3.1 can also be transformed into this form, i.e., replacing gas pressure p and saturated vapor pressure p s with gas density ρ g and adsorbed gas density ρ ad , respectively, and employing (1 − ρ g /ρ ad ) to correct for the true adsorbed gas amount. For example, the Langmuir equation can be transformed into the following form: where b r is a constant similar to the Langmuir constant, which has the relationship ρ L = 1/b r . Langmuir density ρ L refers to the gas density at which the adsorption amount is half of the maximum. The transformation of other models occurs in the same way. A previous study [82] pointed out that gas density is the most meaningful variable in high-pressure sorption areas, so it is recommended that they be used in high-pressure studies instead of pressure.

Adsorption/Absorption Models
Reucroft [83] reported that CO 2 dissolved in coal and caused it to swell in addition to being adsorbed on the coal surface. When studying gas sorption behaviors on polymers, Sato et al. [84] found that gas can not only adsorb on a solid surface, but also can be absorbed into the interior of solid material. Larsen [85] proposed a similar two-component sorption on coal samples. Adsorbed gas and dissolved gas exist on the kerogen surface and in the kerogen interior, respectively, which restrict and connect with each other in shale gas reservoirs [38,41]. To model the two types of sorption (adsorption and absorption), different methods were proposed in previous studies.
The first category is the hybrid type, namely the superposition of gas adsorption law and gas absorption law. Gas adsorption can be described by the abovementioned adsorption models, such as the monolayer adsorption models or the micropore filling models, while gas absorption is described by Henry's law. Meanwhile, supercritical sorption characteristics need to be considered in high-pressure and high-temperature sorption studies. Here, we take the Langmuir adsorption and D-A models as examples to introduce a hybrid method, and other adsorption models, such as the Freundlich, Toth, and D-R models, can be handled by the same procedure.
(1) Langmuir/Henry combination: In this model, gas adsorption is modeled by the Langmuir equation, and gas absorption is described by a term proportional to pressure, following Henry's law. Here, the subcritical adsorption and absorption are described in terms of pressure, as in Equation (28), while supercritical adsorption and absorption are in terms of gas density [81] as in Equation (29): (2) Volume filling/Henry combination: Gas adsorption is described by the D-R or D-A model, and gas absorption is described by Henry's law. For subcritical adsorption and absorption, this can be described in terms of pressure: For supercritical adsorption and absorption description in terms of gas density [81], this is: If gas pressure is adopted in a supercritical adsorption model [86], then the virtual saturated vapor pressure concept introduced in Section 3.1.3 needs to be employed, i.e.: (3) Swelling contribution: Dissolved gas usually swells the solid materials after absorption [83,85]. If the swelling contribution was equal to the condensed gas volume in adsorbents, Equations (28)-(32) need to be modified, because the swelling occupies space that would otherwise be taken up by gases [81]. The term (1 − ρ g /ρ ad ) needs to be multiplied by the absorption term. Taking the supercritical volume filling/Henry combination as an example, the model considering swelling contribution can be expressed as follows: Sakurovs et al. [81] compared the calculation results from the Langmuir model of Equation (27), the D-R model of Equation (26), the Langmuir/Henry combination model of Equation (29), and the D-R/Henry combination model of Equation (31) using CO 2 and CH 4 adsorption data, as shown in Figure 24a,b. The fitting effects of the Langmuir and the D-R models are both poor, while the Langmuir/Henry or D-R/Henry combination model has a much better effect. This means that the added term (kp) in adsorption models improves the fitting effect and reduces the calculated surface sorption capacity. The improvement of this term is greater in the D-R model than in the Langmuir model, since the D-R/Henry combination fits the measured sorption data better than the Langmuir/Henry combination. This suggest that the gas sorption mechanism in coal is more likely the volume filling, rather than monolayer coverage. (4) Bi-Langmuir adsorption model: Assuming absorption and swelling are related, Pini et al. [82] applied the Bi-Langmuir model [87] to describe the combination of adsorption and absorption, where linear superposition is adopted for each Langmuir adsorption term, i.e.: where the first term on the right side of the equation is the adsorption term, and the second term is the absorption term.
Assuming the Langmuir equilibrium constants for adsorption and absorption are equal [82], i.e., b ad = b ab , the excess adsorption amount can be expressed as follows: where G t is the sum of maximum adsorption amount and maximum absorption amount. This Bi-Langmuir model was also compared to the experimental data and the D-R/Henry combination model, as shown in Figure 26. Both models fitted the excess sorption data well, but the fitted curve for absolute sorption data from the D-R/Henry combination model is much higher than the experimental data, while the Bi-Langmuir model fitted the absolute sorption data excellently. This is caused by the neglect of the swelling effect and the assumption of unlimited sorption capacity in the D-R/Henry combination model. Therefore, it should be seen as an empirical approach to describe excess sorption isotherms, and cannot be used for gas storage capacity estimation. Contrarily, the Bi-Langmuir model has a solid physical basis from experimental observations of a saturation-limited equilibrium between the gas phase and the condensed phase [82]. From this point, the Bi-Langmuir model is more suitable for adsorption and absorption modeling of shale gas.

Bulk Gas Properties
Gas Density (1) Calculated by EOS: Generally, the bulk phase gas density can be calculated by the real gas EOS, as mentioned in Section 2.6.1. After calculation, we can also obtain the relationship between gas density and pressure by the nonlinear fitting technique, which expresses density in terms of pressure: where c 0 , c 1 , c 2 , and c 3 are fitting parameters.
As we can see from Figure 27a, free gas density decreases with increasing temperature when 0.1 MPa < p < 30 MPa and increases with increasing pressure when 273.13 K < T < 373.13 K. Compared to temperature, the influence of pressure on free gas density is more obvious. Since the change in pressure is more marked than that of temperature during shale gas reservoir development, more attention needs to be paid to the change in free gas density with pressure. The influence of temperature on free gas density is the most severe at a pressure of around 15 MPa, while it is weaker at lower or higher temperatures.

Bulk Gas Properties
Gas Density (1) Calculated by EOS: Generally, the bulk phase gas density can be calculated by the real gas EOS, as mentioned in Section 2.6.1. After calculation, we can also obtain the relationship between gas density and pressure by the nonlinear fitting technique, which expresses density in terms of pressure: where c0, c1, c2, and c3 are fitting parameters.
As we can see from Figure 27a, free gas density decreases with increasing temperature when 0.1 MPa < p < 30 MPa and increases with increasing pressure when 273.13 K < T < 373.13 K. Compared to temperature, the influence of pressure on free gas density is more obvious. Since the change in pressure is more marked than that of temperature during shale gas reservoir development, more attention needs to be paid to the change in free gas density with pressure. The influence of temperature on free gas density is the most severe at a pressure of around 15 MPa, while it is weaker at lower or higher temperatures. (2) Measured by experiments: Bulk phase gas density can also by measured by experiments. Since analytical modeling is the main method introduced in this article, experimental measurement apparatus and procedures are not introduced.

Gas Viscosity
Accurate determination of natural gas viscosity plays a key role in its management as it is one of the most important parameters in calculations. It is a pressure-and temperature-dependent parameter that can be calculated by different empirical methods. In this section, different calculation formulations  Gas Viscosity Accurate determination of natural gas viscosity plays a key role in its management as it is one of the most important parameters in calculations. It is a pressure-and temperature-dependent parameter that can be calculated by different empirical methods. In this section, different calculation formulations are compared. First, we will introduce some density-based models, where the calculation accuracy of viscosity is based on the prediction of gas density. We note that the following formula in the original work may use different units, and we have transformed all parameters into the SI unit for convenience, i.e., viscosity in Pa s, density in kg/m 3 , molar mass in kg/mol, and temperature in K. With different units, the coefficients in the equation have different values.
(1) Lee method [88,89] Lee et al. reported a viscosity calculation formula of light hydrocarbons based on accurate density data of pure or mixed gas components. Using a linear molecular weight mixing rule, gases' viscosity is expressed as a function of the molecular weight and the gas density, which is: where: where µ is the viscosity, Pa s; M is the molar mass, mol/kg; T is the temperature, K; and K, X, and Y are intermediate parameters.
The results of Equation (38) provided satisfactory fitting capability to the data on methane, ethane, propane, and n-butane simultaneously, with a standard deviation of 1.89% [88,89]. The coefficients in Equations (39)- (41) were determined from the experimental data for pressures ranging from 0.69 MPa to 55.16 MPa and temperatures ranging from 310.93 K to 444.26 K.
(2) Improved Lee method [90,91] If experimental data on gas density are not available, it can be predicted by an empirical method. In this situation, the coefficients of K, X, Y can be obtained as follows: The calculated viscosity of this method agrees with parts of the published viscosity data within 2% at low pressure and within 4% at high pressure when the specific gravity of the gas is smaller than 1.0. The method is less accurate for gases of higher specific gravities, usually giving lower estimates by up to 20% for retrograde gases with specific gravities over 1.5. Different from the original study, we take SI units here in Equations (43)- (45), i.e., µ in Pa·s, ρ g in kg/m 3 , M in kg/mol and T in K.
(3) Londono method [92] Energies 2020, 13 Another polynomial gas viscosity model was proposed in the study of Londono et al. [92] based on nonlinear regression techniques, which can be expressed as: The values of the parameters in Equations (50)-(57) are given in Table 5. Sutton correlated the effect of intermolecular forces as a function of apparent molar weight, pseudocritical pressure, and pseudocritical temperature and used the following coefficients: (6) Abooali-Khamehchi (A-K) method [95] Abooali and Khamehchi developed a natural gas viscosity calculation method covering 1938 data points with the following expression: The above density-based models provide reliable ways to calculate natural gas viscosity [96]. However, the accuracy of the gas density calculation in the models depends on the prediction of gas deviation factor Z at elevated pressure and temperature [97]. To overcome this problem, correlations based on the corresponding states theory were also established.

(7) Heidaryan-Moghadasi-Salarabadi (H-M-S) method [98]
Based on the falling body viscometer experiment, a correlation for methane gas viscosity was presented for temperatures up to 400 K and pressures up to 140 MPa, with an average absolute percent relative error of 0.794: (9) Heidaryan-Esmaeilzadeh-Moghadasi (H-E-M) method [96] In this method, gas viscosity is expressed as follows: where the values of coefficients A e1 -A e10 can be seen in Table 7. A generalized correlation method was proposed based on 3231 data points of 29 multicomponent mixtures at wide ranges of pressure (0.1-137.8 MPa), temperature (241-473 K), and specific gravity (0.573-1.337) by Jarrahian and Heidaryan [97], which can be expressed as follows: where µ 1atm is the gas viscosity at 0.1 MPa; p pr and T pr are pseudoreduced pressure and pseudoreduced temperature; and J 1 , J 2 , and J 3 are fitting parameters with the values J 1 = 7.86338004624174, J 2 = −9.00157084101445 × 10 −6 , and J 3 = 0.278138950019508.
(11) Izadmehr method [100] Two models were developed for pure natural gas and impure natural gas viscosity calculation based on genetic programming techniques, covering 6484 data points and suitable for temperatures ranging from 109.6 K to 600 K, pressures ranging from 0.01 MPa to 199.95 MPa, and gas specific gravity ranging from 0.553 to 1.5741.
For the prediction of sour or sweet natural gases, the following formula can be employed: To improve the precision of pure natural gas viscosity prediction, the viscosity is calculated as follows: where the coefficients of a iz to f iz are shown in Table 8.  (70) and (71). ( As can be seen, gas viscosity at low-pressure conditions, namely 1 atmospheric pressure, is required in some models. Dempsey used graphical correlations and gave the following expression [96]:

Coefficient Sour/Sweet Natural Gas Equation (70) Sweet Natural Gas Equation
Standing improved the procedure of Dempsey for atmospheric viscosity calculation; the result is known as the Dempsey-Standing method [96]: Meanwhile, the coefficient K in the Sutton method above is also equivalent to the low-pressure gas viscosity, the value of which is 10 4 cp in the original work [93]. Figure 28a shows the relationship between gas viscosity at low pressure (p = 0.1 MPa). The Dempsey method [96] leads to much higher values than the other three models, while the Londono method [92] values are slightly higher than those from the Standing method [96] and Sutton method [93]. The Standing method [96] and Sutton method [93] produce basically the same viscosity results. Figure 28b displays gas viscosity variance with temperature at different pressures. For pressure lower than 14 MPa, the gas viscosity increases as temperature increases, which is the opposite situation to that of liquids. With temperature increasing, the gas characteristics become more similar to those of liquids, and gas viscosity decreases as temperature increases.
A correlation method for this calculation was proposed by Londono et al. [92] as follows: Meanwhile, the coefficient K in the Sutton method above is also equivalent to the low-pressure gas viscosity, the value of which is 10 4 cp in the original work [93]. Figure 28a shows the relationship between gas viscosity at low pressure (p = 0.1 MPa). The Dempsey method [96] leads to much higher values than the other three models, while the Londono method [92] values are slightly higher than those from the Standing method [96] and Sutton method [93]. The Standing method [96] and Sutton method [93] produce basically the same viscosity results. Figure 28b displays gas viscosity variance with temperature at different pressures. For pressure lower than 14 MPa, the gas viscosity increases as temperature increases, which is the opposite situation to that of liquids. With temperature increasing, the gas characteristics become more similar to those of liquids, and gas viscosity decreases as temperature increases. In Figure 29, the gas viscosity from different models is compared with the experimental data for pressures ranging from 0.1 MPa to 3.3 MPa at a temperature of 293.15 K. The experimental data were chosen from the study of Hurly et al. [101]. As we can see, the calculation results vary: some match the experimental data well, while others deviate from the experimental data significantly, such as the H-J [94] and H-E-M methods [96]. To quantitatively evaluate the different models, the relative deviation and average absolute relative deviation (AARD), in terms of the experimental data, are given in Figures  30 and 31  In Figure 29, the gas viscosity from different models is compared with the experimental data for pressures ranging from 0.1 MPa to 3.3 MPa at a temperature of 293.15 K. The experimental data were chosen from the study of Hurly et al. [101]. As we can see, the calculation results vary: some match the experimental data well, while others deviate from the experimental data significantly, such as the H-J [94] and H-E-M methods [96]. To quantitatively evaluate the different models, the relative deviation and average absolute relative deviation (AARD), in terms of the experimental data, are given in Figures 30 and 31, respectively. Figure 30 illustrates that the H-J [94] and H-E-M methods [94] have poor performance for gas viscosity prediction for pressures ranging from 0.1 MPa to 3.3 MPa, while the Izadmehr method [100] performs badly for both sour and sweet gas viscosity prediction when the pressure is smaller than 1.4 MPa. This is because this method is obtained based on the regression of different kinds of gases, while the experimental data are for the gas viscosity of pure methane. Although the H-J method [97] predicts gas viscosity accurately when p < 1 MPa, its deviation from the experimental data becomes more and more obvious as pressure increases. The Lee method, Improved Lee method, Londono method, Sutton method, A-K method, Sanjari method, and Izadmehr sweet gas method can predict gas viscosity for 0.1 MPa < p < 3.3 MPa with a relative deviation smaller than 0.3%. Among these methods, the Izadmehr sweet gas, Sanjari, A-K, and Sutton methods perform the best, with AARD equaling to 0.57%, 0.68%, 0.88%, and 0.95%, respectively, as shown in Figure 31.          Figure 32. The results of the H-M-S method vary significantly as the temperature changes. The discrepancies between the different models become more prominent as pressure increases. However, the discrepancies become less obvious as temperature increases.
(a) T = 273. 13 Figure 32. The results of the H-M-S method vary significantly as the temperature changes. The discrepancies between the different models become more prominent as pressure increases. However, the discrepancies become less obvious as temperature increases.  Figure 32. The results of the H-M-S method vary significantly as the temperature changes. The discrepancies between the different models become more prominent as pressure increases. However, the discrepancies become less obvious as temperature increases.
(a) T = 273. 13   According to the relationship between absolute adsorption and excess adsorption, the isotherm linearization method can also be used to obtain adsorbed phase gas density. Since it is not recommended to calculate the virtual saturated vapor pressure, using this method to obtain adsorbed gas density will not be introduced here. One of the conclusions from this paper [77] is that the value of adsorbed gas density is between the critical density (162.66 kg/m 3 , at 4.59 MPa and 190.53 K) and normal boiling density of liquid phase (422.36 kg/m 3 , at 0.101 MPa and 111.67 K), which can be confirmed from the calculation results of different models, as can be seen in Figure 33a. Note that the data of the Ono-Kondo method and ZGR EOS are collected from a previous study of Sudibandriyo et al. [104].
Energies 2020, 13, 5427 43 of 51 recommended to calculate the virtual saturated vapor pressure, using this method to obtain adsorbed gas density will not be introduced here. One of the conclusions from this paper [77] is that the value of adsorbed gas density is between the critical density (162. 66  The adsorbed gas density is independent of temperature in all methods except the Ozawa model [75] of Equation (82), which assumed adsorbed gas as a superheated liquid. It seems that the Ozawa model is more suitable to predict adsorbed gas density since superheated liquid can expand with increasing temperature. However, the excess adsorption, which is calculated by Equation (25) Figure 33. Comparison of adsorbed gas density calculated from different models (a) and the change of (1 − ρ g /ρ a ) with increasing pressure (b).
The adsorbed gas density is independent of temperature in all methods except the Ozawa model [75] of Equation (82), which assumed adsorbed gas as a superheated liquid. It seems that the Ozawa model is more suitable to predict adsorbed gas density since superheated liquid can expand with increasing temperature. However, the excess adsorption, which is calculated by Equation (25), may have a negative value, as can be inferred from Figure 33b. This is not physically right. As we mentioned, ρ b and T b of gases are functions of pressure for superheated liquid. These pressure-dependent characteristics are ignored in the original work [75], which cannot be neglected in a sorption study of shale gas reservoirs. Therefore, corresponding research needs to be done on this subject to develop a more practical model for shale gas sorption study.

Experimental Method
The experimental surface sorption amount can be a starting point for the calculation of absolute adsorption and adsorbed gas density. From Equations (24) and (25), we can obtain another expression of excess sorption: For high-pressure adsorption, bulk gas is more compressible than adsorbed gas and the excess sorption isotherm is mainly influenced by bulk gas EOS, where the adsorbed gas volume can be seen as a constant with increasing pressure. From this standpoint, the line segment after the inflection point on plot of excess sorption amount versus bulk gas density can be used for calculating the volume and density of adsorbed gas, as we can see from Figure 34. The density and volume of adsorbed gas after the maximum sorption amount can be speculated by the slope and intercept of line, according to Equation (83), where the absolute value of the slope is adsorbed gas volume and the intersection point of linear segment and x-axis is adsorbed gas density [80,104]. The calculated adsorbed gas density of methane is between 413.78 kg/m 3 and 433.41 kg/m 3 , which is in the range of liquid density at its boiling point ρ b = 422.36 kg/m 3 . The calculated adsorbed gas volume for methane is between 0.423 cm 3 /g and 0.467 cm 3 /g for 303 K < T < 333 K. Note that, due to data error and the difference in data processing methods, there is a slight variation in the results between our calculation (see in Figure 35) and that of Moellmer et al. [80]. have a negative value, as can be inferred from Figure 33b. This is not physically right. As we mentioned, ρb and Tb of gases are functions of pressure for superheated liquid. These pressure-dependent characteristics are ignored in the original work [75], which cannot be neglected in a sorption study of shale gas reservoirs. Therefore, corresponding research needs to be done on this subject to develop a more practical model for shale gas sorption study.

Experimental Method
The experimental surface sorption amount can be a starting point for the calculation of absolute adsorption and adsorbed gas density. From Equations (24) and (25), we can obtain another expression of excess sorption: For high-pressure adsorption, bulk gas is more compressible than adsorbed gas and the excess sorption isotherm is mainly influenced by bulk gas EOS, where the adsorbed gas volume can be seen as a constant with increasing pressure. From this standpoint, the line segment after the inflection point on plot of excess sorption amount versus bulk gas density can be used for calculating the volume and density of adsorbed gas, as we can see from Figure 34. The density and volume of adsorbed gas after the maximum sorption amount can be speculated by the slope and intercept of line, according to Equation (83), where the absolute value of the slope is adsorbed gas volume and the intersection point of linear segment and x-axis is adsorbed gas density [80,104]. The calculated adsorbed gas density of methane is between 413.78 kg/m 3 and 433.41 kg/m 3 , which is in the range of liquid density at its boiling point ρb = 422.36 kg/m 3 . The calculated adsorbed gas volume for methane is between 0.423 cm 3 /g and 0.467 cm 3 /g for 303 K < T < 333 K. Note that, due to data error and the difference in data processing methods, there is a slight variation in the results between our calculation (see in Figure 35) and that of Moellmer et al. [80].  Excess sorption amount, mg/g Free gas density,kg/m 3 T=333K T=318K T=303K Figure 34. Excess sorption amount of methane in three different temperatures and the line segment after inflection point of maximum [80].  Figure 35. The change in adsorbed gas density (a) and adsorbed gas volume (b) with temperature, calculated by an experimental method.
As we can see from Figure 35a, the adsorbed gas density of methane decreases as temperature increases, showing a similar tendency to the Ozawa prediction [75], but with a higher value. This may provide a way to modify the Ozawa model of Equation (82), as we have mentioned that constant values of ρb and Tb for gases as well as the coefficient −0.0025 may not be suitable for methane sorption studies at high pressure. The calculated adsorbed gas volume of methane also shows a decreasing tendency with increasing temperature, as shown in Figure 35b, which corresponds to a decreasing amount adsorbed at increasing adsorption temperatures [80].

Concluding Remarks
Due to the massive development of shale gas reservoirs in recent years, the understanding of shale formation characteristics and shale gas storage forms has become a research hotspot. Our increased understanding has helped to extend the application of the classical seepage theory to new fields, where nanoscale flow spaces, gas sorption behaviors, and real gas properties are taken into account. However, much disagreement and confusion on this subject still exist, so it requires comprehensive investigation in the future.
This review includes a summary, discussion, and comparison of shale formation characteristics, shale gas occurrence types, and property calculation methods of adsorbed gas and free gas, providing fundamental support to the deep understanding of shale gas reservoirs. (1) The typical mineral composition and organic geochemical characteristics, as well as pore size distribution, of shale Adsorbed gas volume, cm 3 /g Temperature, K Figure 35. The change in adsorbed gas density (a) and adsorbed gas volume (b) with temperature, calculated by an experimental method.
As we can see from Figure 35a, the adsorbed gas density of methane decreases as temperature increases, showing a similar tendency to the Ozawa prediction [75], but with a higher value. This may provide a way to modify the Ozawa model of Equation (82), as we have mentioned that constant values of ρ b and T b for gases as well as the coefficient −0.0025 may not be suitable for methane sorption studies at high pressure. The calculated adsorbed gas volume of methane also shows a decreasing tendency with increasing temperature, as shown in Figure 35b, which corresponds to a decreasing amount adsorbed at increasing adsorption temperatures [80].

Concluding Remarks
Due to the massive development of shale gas reservoirs in recent years, the understanding of shale formation characteristics and shale gas storage forms has become a research hotspot. Our increased understanding has helped to extend the application of the classical seepage theory to new fields, where nanoscale flow spaces, gas sorption behaviors, and real gas properties are taken into account. However, much disagreement and confusion on this subject still exist, so it requires comprehensive investigation in the future.
This review includes a summary, discussion, and comparison of shale formation characteristics, shale gas occurrence types, and property calculation methods of adsorbed gas and free gas, providing fundamental support to the deep understanding of shale gas reservoirs. (1) The typical mineral Energies 2020, 13, 5427 45 of 50 composition and organic geochemical characteristics, as well as pore size distribution, of shale formations are given based on our measurements and analysis of samples from the Longmaxi Formation. We found that mesopores are the mainly developed pore types in the Longmaxi Shale Formation, with the gas-filled porosity of shale samples ranging from 1% to 7.5%, and permeability usually smaller than 0.1 mD. Three shale gas types are usually classified, free gas, adsorbed gas, and dissolved gas, where adsorbed gas and dissolved gas are often considered as one type due to the equilibrium state between them. Therefore, it is claimed that there are two gas types, free gas and adsorbed gas, in shale gas plays in some research. (2) Methane in shale gas reservoirs is in a supercritical state, so Gibbs excess sorption models and supercritical state sorption models are employed to capture the gas adsorption and desorption behaviors. Meanwhile, different models considering not only the gas adsorption but also the absorption are introduced. Great discrepancies could occur if the supercritical state and Gibbs excess adsorption characteristics are ignored. Different mechanisms of adsorption in micropores and macropores may explain the hysteresis between adsorption and desorption, rather than capillary condensation. (3) Different methods of calculating gas properties, such gas free gas density, free gas viscosity, and adsorbed gas density considering high-pressure and high-temperature conditions, are summarized, with recommended approaches given after the comparison. From our review, we can see that the geological characteristics of shale formations are quite different from those of conventional ones, and need further assessment using a high-resolution apparatus. Gas adsorption mechanisms are still not clear, although numerous models have been developed to account for this phenomenon. For example, current supercritical adsorption models based on the potential theory are modifications of a previous adsorption theory, which are dependent on empirical parameters and lack of a universal theoretical basis. Applicable multicomponent gas adsorption models considering high-pressure and high-temperature conditions and pore size distribution are in demand to describe gas behavior in shale gas reservoirs. They are a rough way to clarify adsorbed gas and dissolved as one type, since dissolved gas may play an important role in shale gas production. Therefore, current adsorption and absorption models should be improved, based on the practical relationship between adsorbed gas and dissolved gas under in situ conditions.