Engineering, Technology & Applied Science Research Vol. 7, No. 4, 2017, 1843-1849 1843 www.etasr.com Rajabi et al.: A Numerical Study of Gas Injection and Caprock Leakage from Yort-e-Shah Aquifer in Iran A Numerical Study of Gas Injection and Caprock Leakage from Yort-e-Shah Aquifer in Iran Maysam Rajabi Faculty of Mining and Metallurgical Engineering Amirkabir University of Technology Tehran, Iran Hossein Salari Rad Faculty of Mining and Metallurgical Engineering Amirkabir University of Technology Tehran, Iran Mohsen S. Masoudian Nottingham Centre for Geomechanics, University of Nottingham University Park, United Kingdom Abstract—In order to mitigate the adverse effects of global warming due to anthropogenic CO2 emission into the atmosphere, geological sequestration of CO2 into subsurface formations has been investigated by many studies over the last decade. However, selection of formations and sites for any field application is still open to debate. The most important properties of a formation suitable for carbon sequestration are those which impact the fluid flow processes. The injection or extraction of gas can change the pore pressure within the reservoir, which in turn results in redistribution of the stress field. These events may cause considerably leakage of the fluid into the surrounding geological formations or ground surface. The main objective of this paper is to evaluate the potential of Yort-e-Shah aquifer for CO2 storage, through a series of analyses with a simplified numerical model. The numerical results suggest that the optimum injection pressure in Yort-e-Shah aquifer is about 15.51 MPa with a safety factor of about 1.7. The results of the fluid pressure and gas plume expansion are presented. Also, an analysis was carried out for a case with leak through cap rock. When there is no leak, the pressure within the aquifer is stable, while on the other hand, the pressure in case of leakage is slightly smaller. In case of leakage, the pressure is lowest in the middle of the reservoir, mainly because the nodes at the middle of the aquifer are influenced by all the leakage points, while around the wellbore or near the end of gas plume, are affected less due to their longer distance to leakage points. Keywords-Numerical model; Fluid flow; Co2 sequestration; Optimum gas pressure; leakage through cap rock I. INTRODUCTION Global warming has a direct effect on climate patterns and this includes the increasingly frequent occurrence of extreme weather events such as heat waves, floods, storms, droughts and bushfires [1,2]. The increase in the global surface temperature from 1956 to 2005 is 0.13 oC per decade and eleven of the twelve years between 1995 and 2006 rank among the twelve warmest years since 1850 [3, 4]. Scientists have suggested several options through which the amount of greenhouse gas emissions into the atmosphere can be reduced. Among these, geological sequestration of carbon dioxide has received a lot of attentions over the last decade [5, 6], and different formations may be suitable for CO2 storage. The estimated capacity of CO2 storage in geological formations are listed in Table I, where it can be seen that saline aquifers offer the highest potential storage capacity. The feasibility of storage in aquifers has been extensively studied [7] and it has exhibited some technological and economic attractiveness. The main parameters in selecting an aquifer for CO2 storages, are its size, porosity, permeability and depth [8, 9]. In addition, an overlying formation is required to cover the aquifer which provides an impermeable seal to prevent the leakage of gas into the atmosphere. The injected CO2 stored in aquifer by a combination of three mechanisms, namely immobilization in traps, dissolution in water, and geochemical interaction with reservoir rocks. In order to assess the suitability of the reservoir for gas storage purposes, a great number of theoretical, numerical, and analytical approaches have been employed [10- 14] and strategies for site selection have been proposed [15]. Since the first CO2 injection into saline aquifer (in the 90s), many progresses have been made and many pilot and filed scale projects have been planned and/or are being undertaken. TABLE I. GLOBAL STORAGE CAPACITY FOR GEOLOGICAL SEQUESTRATION OPTIONS Reservoir type Storage capacity (Gt CO2) Lower estimate Upper estimate Oil and gas fields 675a 900a Un- mineable coal seams (ECBM) 3-15 200 Deep saline formation 1000 Uncertain,but possibly 10000 a These numbers would increase by 25% if undiscovered oil and gas were included in this assessment In Iran, dating back at the 70s, different areas were studied to shortlist suitable locations. Based on the results of the primary studies, the Yort-e-Shah Region in Varamin was chosen for further studies by the National Iranian Oil Company [16]. Therefore, Yort-e-Shah aquifer can be alternatively considered for gas storage. This paper provides the results of numerical simulations where the potential of Yort-e-Shah aquifer for carbon sequestration is investigated. In order to achieve this, a simplified numerical model was employed where the radial flow is considered around a single vertical wellbore. Through a series of simulations the optimum injection pressure and the storage capacity of Yort-e-Shah aquifer is determined. Also, an analysis is carried out for a case Engineering, Technology & Applied Science Research Vol. 7, No. 4, 2017, 1843-1849 1844 www.etasr.com Rajabi et al.: A Numerical Study of Gas Injection and Caprock Leakage from Yort-e-Shah Aquifer in Iran with leaking through cap rock. It should be noted that the model taking into account the numerical study of CO2 injection and leakage through caprock, has been modified and verified compared with the previous study [17, 18]. The results of this study, however, are advantageous to many projects involving gas injection into subsurface formations, enhanced oil and gas recovery. II. GEOLOGICAL STUDIES OF YORT-E-SHAH AQUIFER IN VARAMIN Yort-e-Shah aquifer is in the form of an anticline and is located 70 kilometers southeast of Tehran and 40 kilometers south of Varamin. Yort-e-shah anticline was formed in 4 stages by extension during middle Alpine movements, then stagnation and reactivated extension by late middle Alpine movements, and finally compressions by late Alpine movements [16]. Detailed geological studies of this aquifer, started in 1997, for the purpose of natural gas storage. The primary phase of the project carried out by establishing descriptive wellbores number 2 and 3 by KBB. Then, Sofregaz was responsible for drilling wells number 4 and 5, and also performed supplementary studies. From the geological point of view, there are three formations in the region as listed in Table II. The Upper Red formation is a sequence of sandstone, clay stone and anhydrite, below which Qum formation contains sequences of anhydrite, gypsum and limestone. Below Qum formation, there is a Volcanic formation. Qum formation has shown promising properties for storage of gas, while the Upper Red formation has shown high potential to prevent gas leakage and is considered as a suitable cap rock, given its very low permeability and porosity measured from laboratory analysis, as listed in Table III [16]. III. MODELING In order to model and evaluate the capacity of storage in Yort-e-Shah aquifer, a simplified numerical model is employed. The model considers a vertical wellbore in the middle of the aquifer, through which gas can be injected. The gas plume, which is the only phase and the component, can flow in the radial and angular direction through the uniform thickness of the aquifer, while the cap rock provides a seal to the storage site. This ideal conditions are depicted in Figure 1. TABLE II. PROPERTIES OF YORT-E-SHAH GEOLOGICAL FORMATIONS Upper depth in well No. 4 (m) Upper depth in well No. 3 (m) Upper depth in well No. 2 (m) Geological sequence Formation name 33 20.50 16 Sequence of sandstone, clay stone and anhydrite Upper Red Formation 948 1002.50 575.10 Sequence of anhydrite, limestone and gypsum Qum Formation 1273 1396.8 804.9 Porphyry andesite Volcanic Formation TABLE III. LABORATORY ANALYSIS OF CAPROCK Sample number Depth (m) Brine Permeability (m2) Gas drive pressure (MPa) Permeability(m2) 14V 571.45 7.007∙10-21 3.445 No flow 5.516 No flow 6.895 No flow 11.032 No flow 13.790 No flow 17.926 No flow 26.890 1.974∙10-23 20V 585.8 3.948∙10-21 3.445 No flow 5.516 No flow 6.895 No flow 11.032 No flow 13.790 No flow 17.926 No flow 26.890 No flow 29V 618.75 4.935∙10-20 3.445 No flow 5.516 No flow 6.895 No flow 11.032 No flow 13.790 No flow 17.926 No flow 26.890 <9.869∙10-23 39V 635.70 1.875∙10-20 3.445 No flow 5.516 No flow 6.895 No flow 11.032 No flow 13.790 No flow 17.926 No flow 26.890 <9.869∙10-23 Fig. 1. Idealized schematic illustration of gas storage in the aquifer. The outer periphery of the aquifer is impervious to gas flow, as is the underlying rock formation. There is an axial symmetry around the single vertical injection wellbore. The well is subject to a variable pressure which was known as function of time. The viscosity of the gas is assumed to be constant, while its density slightly increases with pressure, i.e. slightly compressible gas. It is also assumed that the flow of gas is an isothermal process, an assumption which is made in many of reservoir simulation works. In order to assess the conditions under which the imperfection of the cap rock seal leads to gas leakage from the reservoir, a leakage rate has also been employed. The leakage rate is assumed to only depend on the local pressure of the gas plume and its spatial position. To describe the mathematical foundation and governing equations for gas flow, we start with Darcy’s law as below: r k p v r     (1) Engineering, Technology & Applied Science Research Vol. 7, No. 4, 2017, 1843-1849 1845 www.etasr.com Rajabi et al.: A Numerical Study of Gas Injection and Caprock Leakage from Yort-e-Shah Aquifer in Iran k p v       (2) Where p is local pressure in the gas plume, k is the intrinsic permeability of the aquifer, μ is gas viscosity, vr is superficial radial velocity component, νθ is superficial angular velocity component. The equation of continuity, taking into account radial and angular flows, accumulation due to density change and leakage through the cap rock, gives: t mv r rv rr              )( 1 ),( 1 (3) where ρ is the gas density, φis the porosity of the aquifer and t is the injection time. Note that m is the leakage parameter that defines the mass flow rate of gas per unit surface area of cap rock per unit depth of the aquifer. In other words, the leakage rate has been added as a sink term into the flow equation. Following [19], parameter m is defined as a function of gas pressure as below: ll ppppkm  ),( * (4) lppm  ,0 (5) Where pl is the threshold pressure beyond which the leakage through cap rock occurs, and k* is a coefficient that defines how the leakage rate depends on the gas plume pressure. In other words, the leakage rate is dined as a piecewise function whose value is zero when the gas pressure is below the leakage threshold, and it is linearly related to the differential pressure across leakage, when the gas pressure is beyond the leakage threshold. The gas is considered to be slightly compressible and hence follows this simple equation: cp (6) Where c is the rock compressibility. Substituting Darcy’s law into the equation of continuity yield the basic flow equation as below: t p kkc mp r p r p r p r p r p r p p                              2 22 2 2 2 2 2 1 (7) It is convenient to introduce the dimensionless flow equation, using the following non-dimensional variables: 0p p p  (8) 1 0 r tkp T   (9) 2 0 2 1 kcp rm M   (10) 1r r  (11) where, p0 is the initial uniform pressure throughout the aquifer, r1 is the radius of injection well and r2 is radius of gas plume within the aquifer. Equation (7) will then become: T P M PpPPPpP P                            2 22 2 2 2 2 2 1  (12) Equation (12) is a non-linear second-order partial differential equation that mathematically explain the transient diffusion of pressure within the radial and angular positions distance within the aquifer. This equation can be solved using a finite difference scheme. For the purpose of a more convenient finite difference solution, we can change the radial distance variable into the following logarithmic variable:       1 2ln r rR (13) The flow equation then becomes as follows: T P UQ PP P R P R P P                         2 2 22 2 2  (14) Where RMeTRQQ 2),,(   (15) ReRUU 2)(  (16) To simplify this relationship, equation (14) can be solved for one vertical half of the gas plume, i.e. in the semi-circular region. Both the curved boundary at R=ln(r2/r1), i.e. the gas/water interface, and the meridian plane represented by θ=0 and θ=π will be effectively impervious to gas flow. An approximation to the solution of (14) and its associated initial and boundary conditions can be obtained by a finite difference technique. We introduce a series of grid points spaced uniformly by a dimensionless radial increment ΔR and an angular increment Δθ. Subscripts i and j may then be used to denote that grid point having coordinates R=iΔR and θ=jΔθ. Details of this model are presented in [19]. In order to verify the numerical model, its results can be compared with the results of an analytical model. Assuming an incompressible fluid is injected at constant rate, Q, into a horizontal homogenous isotropic reservoir with an infinite extent and initial pressure of P0, through a vertical wellbore at its center, an analytical solution can be found following [6].              fff ff f r c tk r Ei Hk Q PtrP    44 ),( 2 0 (17) Where H is thickness of the aquifer and Ei is the exponential integral function, which can also be defined for any positive y using Ei, as below:     1 1 )()( dyy e yEiyEi y (18) An example was analyzed by both the numerical and analytical solutions, and the pressure predictions exhibit perfect match, as depicted in Figure 2. IV. RESULTS AND DISCUSSION The main problem investigated in this study, is to determine the optimum injection pressure and time needed to fill the aquifer with gas. As such, the potential capacity of gas storage in the aquifer can be assessed. The optimum injection pressure Engineering, Technology & Applied Science Research Vol. 7, No. 4, 2017, 1843-1849 1846 www.etasr.com Rajabi et al.: A Numerical Study of Gas Injection and Caprock Leakage from Yort-e-Shah Aquifer in Iran and time, can be estimated following the concept developed in [19]. In order to do that, the properties of Yort-e-Shah aquifer is used as an input for the numerical model, as listed in Table IV. The case was then simulated for different gas injection pressures. For each case, the injection time required for the gas pressure within the aquifer to nearly reach the injection pressure is estimated. The injection pressure under which the shortest injection time, with suitable safety factor, was achieved considered to be the optimum injection pressure. This approach was employed for wellbores number 2 to 4 in Yort-e- Shah aquifer. The results of these analyses are shown in Table V. TABLE IV. PROPERTIES OF YORT-E-SHAH AQUIFER USED IN NUMERICAL SIMULATION Parameters Symbol Value Units Reference Averaged porosity in wellbore No. 2 φf 0.086 - [16, 20] Averaged porosity in wellbore No. 3 φf 0.069 - [21, 22] Averaged porosity in wellbore No. 4 φf 0.08 - [18, 23] Averaged permeability in well wellbore No. 2 kf 1.59∙10 -15 m2 [16, 20] Averaged permeability in wellbore No. 3 kf 5.05∙10 -15 m2 [21, 22] Averaged permeability in wellbore No. 4 kf 2.44∙10 -15 m2 [18, 23] Gas viscosity μf 1∙10-4 Pa s [24] Reference gas density 0f 955.3 Kg/m3 From real gas equation of state Initial reservoir pressure pi 9.7 MPa [21, 22] Radius of the reservoir r0 500 m [21, 22] Radius of the wellbore rw 0.13 m [21, 22] Reservoir depth in wellbore No. 2 dr 804.9 m [21, 22] Reservoir depth in wellbore No. 3 dr 1396.8 m [21, 22] Reservoir depth in wellbore No. 4 dr 1273 m [18, 23] Reservoir thickness wellbore No. 2 tr 229.8 m [21, 22] Reservoir thickness in wellbore No. 3 tr 394.3 m [21, 22] Reservoir thickness in wellbore No. 4 tr 325 m [18, 23] Numerical time step Δt 0.1 days - TABLE V. RELATIONSHIP BETWEEN GAS PRESSURE AND INJECTION TIME Well No. 4 Injection Time (day) Well No. 3 Injection Time (day) Well No. 2 Injection Time (day) Injection Pressure (MPa) 635 184 1195 13.448 623 165 1170 13.789 610 152 1145 14.478 595 145 1115 15.168 586 141 1090 15.513 612 175 1160 15.858 597 152 1120 16.547 585 145 1095 17.237 579 150 1095 18.616 555 140 1065 19.305 From the data presented in this table, it can be seen that increasing pressure generally leads to shorter injection time, regardless of injection well. To determine the optimum injection pressure, however, it is essential to consider the leakage threshold. The leakage threshold was assumed to be 27 MPa (according to the Table III). Fig. 2. Pressure profiles at different times predicted by the analytical solution and numerical simulation In order to reduce the risk of leakage a safety factor of 1.7 is chosen, based on which, the optimum gas injection pressure is determined to be 15.5 MPa. Assuming no leakage through the cap rock, the simulations were repeated for each depth within the reservoir, considering the corresponding Engineering, Technology & Applied Science Research Vol. 7, No. 4, 2017, 1843-1849 1847 www.etasr.com Rajabi et al.: A Numerical Study of Gas Injection and Caprock Leakage from Yort-e-Shah Aquifer in Iran permeability and porosity and an injection pressure of 15.5 Mpa applied for 330 days. Figures 3-5 depict the typical variations of gas pressure as a function of depth for the storage wells No. 2 to 4. It can be seen that for wellbore No. 4, the maximum gas pressure within the reservoir was achieved at a depth of about 1034 m. Also, it can be seen that in the maximum pressure was achieved at a depth of 623 m and 1344m for wellbore No. 2 and wellbore No. 3, respectively. From the results given in Figures 3-5, no relationship between depth, porosity, permeability and maximum injection pressure can be concluded, as this is mainly due to the variability of the reservoir. To study the effect of leakage through cap rock, on pressure distribution within the reservoir, a case was considered in which a cascade of leaking points contributed to the leakage (Figure 6). In this case, injection was simulated at a depth of 1041.15 m in wellbore No. 4, while the leakage threshold was assumed to be pL=14.8 MPa and the leakage coefficient was considered to be k*=0.0047 kg/cu.m.day.pa. The same simulation was performed without leakage and the results are summarized in Figure 7. The results for the case with leak are also summarized in Figure 8. As expected, when there is no leak, the pressure within the aquifer is stable, while on the other hand, the pressure in case of leakage is slightly smaller. In case of leakage, the pressure is lowest in the middle of the reservoir, mainly because the nodes at the middle of the aquifer are influenced by all the leakage points, while around the wellbore or near the end of gas plume, are affected less due to their longer distance to leakage points. Fig. 3. Gas pressure, porosity and permeability as a function of depth at wellbore No. 2 Fig. 4. Gas pressure, porosity and permeability as a function of depth at wellbore No. 3 Engineering, Technology & Applied Science Research Vol. 7, No. 4, 2017, 1843-1849 1848 www.etasr.com Rajabi et al.: A Numerical Study of Gas Injection and Caprock Leakage from Yort-e-Shah Aquifer in Iran Fig. 5. Gas pressure, porosity and permeability as a function of depth at wellbore No. 4 Fig. 6. Location of leak in the aquifer Fig. 7. Gas pressure within aquifer at depth of 1041.15 m in wellbore No. 4 (no leak) Fig. 8. The gas pressure in depth 1041.15 m of wellbore No. 4(leaky condition) V. CONCLUDING REMARKS To reduce the adverse effects of global warming, geological sequestration has been suggested, which includes capturing and pumping anthropogenic The Yort-e-Shah aquifer near Tehran, was investigated in this study as a potential candidate for gas storage. One of the main practical aspects for gas storage is to determine the optimum gas pressure and injection time. No particular relationship was found to connect the optimum injection pressure with the depth of the injection within the reservoir. The results of this study showed that an optimum pressure of about 15.5 MPa, provides a leakage safety factor of 1.7. A comparative analysis was also conducted to study the effect leakage on pressure distribution within the aquifer. The leakage, as expected, resulted to a lower pressure level within the aquifer. To determine the suitability of this site for gas storage or CO2 sequestration, further studies are needed, including a more rigorous parametric study. In addition, more exploratory wellbores are needed to understand the extent of variability in the reservoir characteristics. In addition, to evaluating the leakage potential through cap rock, more comprehensive field and experimental investigations are required, where the hydro-mechanical properties of the cap rock are examined. The results of this preliminary study provide some useful first-hand estimation of the gas pressure distribution within the reservoir, with or without leakage potential. The model provided in this paper, can also be further Engineering, Technology & Applied Science Research Vol. 7, No. 4, 2017, 1843-1849 1849 www.etasr.com Rajabi et al.: A Numerical Study of Gas Injection and Caprock Leakage from Yort-e-Shah Aquifer in Iran developed to include the geomechanical response of the reservoir to gas injection. REFERENCES [1] J. Hansen, D. Johnson, A. Lacis, S. Lebedeff, P. Lee, D. Rind, G. Russell, “Climate impact of increasing atmospheric carbon dioxide”, Science, Vol. 213, No. 4511, pp. 957-966, 1981 [2] A. Dai, “Drought under global warming: a review”, Wiley Interdisciplinary Reviews: Climate Change, Vol. 2, No.1, pp. 45-65, 2011 [3] D. A. Lashof, D. R. Ahuja, “Relative contributions of greenhouse gas emissions to global warming”, Nature, Vol. 344, No. 6266, pp. 529-531, 1990 [4] R. K. Pachauri, A. Reisinger, Contribution of working groups I, II and III to the fourth assessment report of the Intergovernmental Panel on Climate Change (IPCC), Geneva, 2005 [5] M. S. Masoudian, “Multiphysics of carbon dioxide sequestration in coalbeds: A review with a focus on geomechanical characteristics of coal”, Journal of Rock Mechanics and Geotechnical Engineering, Vol. 8, No. 1, pp. 93-112, 2016 [6] M. S. Masoudian, D. Airey, A. El-Zein, “Mechanical and flow behaviours and their interactions in coalbed geosequestration of CO2”, Geomechanics and Geoengineering: An International Journal, Vol. 8, No. 4, pp. 229-243, 2013 [7] J. K. Eccles, L. Pratson, R. G. Newell, R. B. Jackson, “Physical and Economic Potential of Geological CO2 storage in saline aquifers”, Environmental Science & Technology, Vol. 43, No. 6, pp. 1962-1969, 2009 [8] M. Bentham, G. Kirby, “CO2 storage in saline aquifers”, Oil & Gas Science and Technology – Revue d'IFP Energies nouvelles, Vol. 60, No. 3, pp. 559-567, 2005 [9] IPCC (Intergovernmental Panel on Climate Change), The IPCC special report on carbon dioxide capture and storage, Cambridge: Cambridge University Press, 2005 [10] S. Bachu, J. M. Nordbotten, M. A. Celia, “Evaluation of the spread of acid gas plumes injected in deep saline aquifers in western Canada as an analogue to CO2 injection in continental sedimentary basins”, In Proceedings of 7th ICGGCT, Vol. 1, 2005 [11] K. Pruess, “On CO2 fluid flow and heat transfer behavior in the subsurface following leakage from a geologic storage reservoir”, Environmental Geology, Vol. 54, No. 8, pp. 1677-1686, 2008 [12] S. L. Bryant, S. Lakshminarasimhan, G. A. Pope, “Buoyancy-dominated multiphase flow and its impact on geological sequestration of CO2”, SPE Journal, Vol. 13, No. 4, pp. 447-454, 2006 [13] Z. Zhang, R. K. Agarwal, “Numerical simulation and optimization of CO2 sequestration in saline aquifers for vertical and horizontal well injection”, Computational Geoscience, Vol. 16, No. 4, pp. 891-899, 2012 [14] L. Nghiem, P. Sammon, J. Grabenstetter, H.Ohkuma, “Modeling CO2 storage in aquifers with a fully-coupled geochemical EOS compositional simulator”, SPE/DOE Symposium on Improved Oil Recovery, Tulsa, April 17–21, p. 51–67, 2004 [15] S. Bachu, “Sequestration of CO2 in geological media in response to climate change: road map for site selection using the transform of the geological space into the CO2 phase space”, Energy Conversion and Management, Vol. 43, No. 1, pp. 87-102, 2002 [16] D. Powell, Lab investigation of Yort-e-Shah aquifer: Petrophysical core analysis study of exploration well No. 2, Task Report A10, Vol. 3, KBB Group, Hannover, 1998 [17] M. Rajabi, H. Tavakoli, M. Shafiei, “Prediction of optimum gas pressure and injection time in underground gas storage reservoirs”, Proceedings of 3rd Iranian Mining Engineering Conference, Yazd, January 26–29, p. 72–79, 2010 [18] P. Werner, Drilling program for Yort-e-Shah well No. 4, Task Report A2.1, Sofregaz Group, Paris, 2006 [19] E. Brooks, J. Ellis, G. Grow, E. Hedges, C. Stout, T. Walsh, New Concepts in Underground Storage of Natural Gas, Report PO-50, University of Michigan, Ann Arbor, 1966 [20] K. Wiechart, Lab investigation of Yort-e-Shah aquifer: Petrophysical core analysis study of exploration well No. 2, Task Report A10, Vol. 2, KBB Group, Hannover, 1998 [21] D. Powell, Lab investigation of Yort-e-Shah aquifer: Petrophysical core analysis study of exploration well No. 3, Task Report A10, Vol. 3, KBB Group, Hannover, 1999 [22] K. Wiechart, Lab investigation of Yort-e-Shah aquifer: Petrophysical core analysis study of exploration well No. 3, Task Report A10, Vol. 2, KBB Group, Hannover, 1999. [23] M.R. Esfahani, E. Kazemzadeh, J. Vali, Routine core analysis report for well No. 4 of Yort-e-Shah field, Core Research Department, Exploration & Production Research Division, Research Institute of Petroleum Industry (RIPI), Tehran, 2006 [24] A. Fenghour, W.A. Wakeham, V. Vesovic, “The viscosity of carbon dioxide”, Journal of Physical Chemistry Reference Data, Vol. 27, No. 1, pp. 31–44, 1998