Table of Contents Author Guidelines Submit a Manuscript
Volume 2017, Article ID 6803294, 11 pages
Research Article

Application of Numerical Modelling and Genetic Programming in Hydrocarbon Seepage Prediction and Control for Crude Oil Storage Unlined Rock Caverns

1Mineral Industries Research Center, Shahid Bahonar University of Kerman, Kerman, Iran
2Mining Engineering Department, Shahid Bahonar University of Kerman, Kerman, Iran
3Chemical Engineering Department, Shahid Bahonar University of Kerman, Kerman, Iran

Correspondence should be addressed to Ebrahim Ghotbi Ravandi; moc.oohay@miharbe_ibtohg

Received 23 November 2016; Revised 2 March 2017; Accepted 15 March 2017; Published 14 June 2017

Academic Editor: Micol Todesco

Copyright © 2017 Ebrahim Ghotbi Ravandi et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.


Seepage control is a prerequisite for hydrocarbon storage in unlined rock caverns (URCs) where the seepage of stored products to the surrounding host rock and groundwater can cause serious environmental and financial problems. Practically seepage control is performed by permeability and hydrodynamic control methods. This paper employs numerical modelling and genetic programming (GP) for the purpose of seepage prediction and control method determination for the crude oil storage URCs based on the effective parameters including hydrogeologic characteristic of the rock and physicochemical properties of the hydrocarbons. Several levels for each parameter were considered and all the possible scenarios were modelled numerically for the two-phase mixture model formulation. The corresponding seepage values were evaluated to be used as genetic programming data base to generate representative equations for the hydrocarbon seepage value. The coefficients of determination () and relative percent errors of the proposed equations show their ability in the seepage prediction and permeability or hydrodynamic control method determination and design. The results can be used for crude oil storage URCs worldwide.

1. Introduction

Underground storage of hydrocarbons in unlined rock caverns (URCs) is more secure, safe, and economical than above-ground storage and has several environmental and operational advantages. The main concern associated with URCs is the seepage of stored products, vapor, and VOCs of them to the surrounding host rock and groundwater which can cause serious environmental problems such as groundwater contamination and accumulation of flammable gas near the surface as well as economic and financial losses. Therefore, seepage control is a fundamental prerequisite for URCs where minimum product seepage is required. Practically seepage control from an URC is performed by permeability or hydrodynamic control methods. Permeability control means applying techniques such as grouting or freezing to control and decrease hydrocarbon seepage by maintaining a specified low permeability and sealing of the rock mass. However these techniques are very time-consuming and expensive [1]. By hydrodynamic control, it is meant that there is groundwater in the rock mass with the static head that exceeds the internal storage pressure resulting in positive groundwater gradient towards the cavern to prevent the escape of the stored products [2]. Aberg [3] postulated that no seepage will occur if the water pressure gradient towards the cavern is positive and greater than unity. There is no standard for acceptable seepage value and how much sealing work is required. It depends on the environmental and economic (operational) aspects. As a general rule, 24 m3/24 hr period in a cavern of 100,000 m3 is considered to be an acceptable limit [1]. Liquid hydrocarbons (e.g., crude oil and gasoline) are not stored under pressure and their pressure inside caverns is hydrostatic. Hydrocarbons vapors and gases pressure inside caverns varies as a function of the temperature, oil level, and components, usually from 0.5 to 3 bar (105 Pa) [4]. Typically, 75% and 25% of a cavern space are considered for liquid oils and gases, respectively [2].

Hydrocarbons seepage (especially crude oil) from underground unlined rock caverns to the surrounding saturated porous media is barely investigated in the literature [5, 6]. In this paper prediction equations for the hydrocarbons (crude oil and gas) seepage value from the Iranian URCs in terms of m3/24 hr are represented using genetic programming based on the data gathered by numerical modelling of the hydrocarbon seepage for a variety of conditions using COMSOL and the two-phase mixture model formulation (see Supplementary Material available online at By applying and solving the proposed equations for the seepage values less than the allowable one, prejudgment can be done and the seepage control technique (e.g., permeability or hydrodynamic control) can be selected.

2. Two-Phase Mixture Model

Two-phase mixture model which was first mentioned by Wang and Beckermann [7] uses mixture variables to reduce the number of partial differential equations (PDEs) of classical two-phase fluid flow formulations in porous media. Therefore, it is more convenient to use the appropriate numerical schemes for the two-phase mixture model [7].

In the mixture model, the two phases are regarded as constituents of a binary mixture and the mixture variables such as mixture density and mixture velocity vector are denoted without index. With introduction of the mixture quantities (see [7]) the conservation of mass with the porosity () is defined asConservation of momentum using Darcy velocity is as follows in which the dynamic viscosity and the pressure are also mixture quantities (see [7]):where is the kinetic mixture density, is the intrinsic permeability, is the saturation and the subscripts of and are related to the wetting and nonwetting phase, respectively, is the relative permeability, is the gravitational acceleration vector, and is the depth.

The diffusive mass flux connects the mixture mean velocity with the velocity of the individual phases:Wang and Beckerman [7] introduced the diffusive flux () as follows:where is the kinematic viscosity, is the hindrance function for phase migration and separation, and is the capillary diffusion coefficient as a function of the wetting phase saturation:where is capillary pressure. The transport equation is as follows, where is the fluid content of the wetting phase:Several scientists have tried to derive a functional correlation for the relative permeability and the capillary pressure () as a function of the wetting phase saturation based on the experimental data. Brooks and Corey (1964) developed an empirical correlation utilizing the entry capillary pressure () and the wetting phase saturation that empirically has been found to be appropriate for the drainage process as follows [8, 9]:where is the pore size distribution index. Its value is usually considered to be 2 for the carbonated rocks [10].

Exact solutions for two-phase fluid flow problems in porous media which involve gravity, capillarity, and fluid flow in two or three dimensions (multidimensional flow) are impossible due to inherent nonlinearity and the need to solve for multiple dependent variables along with a variety of unknowns. Solving practical problems requires a suitable numerical method [7]. A lot of authors have used numerical methods and software tools to model single- or two-phase fluid flow in porous and fractured media [1115]. In the mentioned literature the effect of gravity, capillary pressures, and multidimensional flow is usually neglected and not considered simultaneously.

3. Validation of Numerical Modelling of Two-Phase Mixture Model by COMSOL

Neglecting the capillary pressure and gravity effects, the five-spot problem is the standard porous media problem where a square computational domain is initially saturated with the nonwetting phase (oil) and the wetting phase (water) is injected through a well at a lower corner of the domain at a constant rate (or pressure) and displaces the oil. The nonwetting phase is produced at the same rate through a well in the opposite upper corner. In order to evaluate the computational efficiency and accuracy of the mixture model formulation numerical modelling with COMSOL, a verification example of the computational domain with the dimensions 300 m × 300 m for the five-spot problem is given to compare the numerical results with the fully coupled (classical) formulation. Boundary and initial conditions are depicted in Figure 1. Dirichlet pressure boundary conditions are 5 m × 5 m injection and production wells. The other boundaries are impermeable and Neumann no-flow boundary conditions. The nonwetting and wetting phase density and viscosity are considered 1000 Kg/m3 and 0.001 Pa·m, respectively. Intrinsic permeability, porosity, and pore size distribution index are 10−7 m2, 0.2, and 2, respectively. Figure 1 shows the comparative study of the numerical modelling results for the fully coupled [16] and mixture model formulations referring to the time days and time step of 1 day. As it can be seen from Figure 1 the results (wetting phase saturation fronts and contour lines) of the fully coupled and mixture model formulations are in good agreement with each other.

Figure 1: Initial and boundary conditions for the five-spot problem as well as water saturation () fronts and contour lines for the time days and the time step () of 1 day for the fully coupled and mixture model formulation with COMSOL.

4. Genetic Programming

Genetic programming (GP) as an extension of the genetic algorithms (GA) was introduced by Koza [17]. The main difference between genetic programming and genetic algorithms is the representation of the solution. Genetic algorithms create a string of numbers that represent the solution but genetic programming creates computer programs (CPs) in the lisp or scheme computer languages as the solution [18]. GP can be used to find a relationship between input and output data in the form of mathematical expression represented by the functions generated in the training (learning) process. If the error rate reaches a certain threshold, the training can be stopped and the testing (validation) can be applied to verify the effectiveness of the best function.

Genetic programming uses four steps to solve problems [18]:(1)Generate an initial population of random compositions of the functions and terminals of the problem (computer programs).(2)Execute each program in the population and assign it a fitness value according to how well it solves the problem.(3)Create a new population of computer programs by applying the following genetic operations:(i)Copy the best existing programs (reproduction).(ii)Create new computer programs by mutation.(iii)Create new computer programs by crossover.(4)The best computer program that appeared in any generation (the best-so-far solution) is designated as the result of genetic programming.Figure 2 represents the genetic programming flowchart.

Figure 2: Genetic programming flowchart.

5. Methodology

The purpose of this study is to employ numerical modelling and genetic programming to predict the hydrocarbon seepage from the URCs (Iranian URCs in the limestone rocks) based on the allowable seepage value (m3/24 hr) to be able to decide on the seepage control technique selection. To reach the stated goal, two parts of numerical modelling were done for the oil and gas seepage from the URCs. The influencing parameters on the hydrocarbons seepage including hydrogeological properties of the rock mass and physicochemical properties of the hydrocarbons were considered in several levels and using full factorial design all the hypothetical cases were modelled using the finite element based commercial software COMSOL version 5.1 and the two-phase mixture model formulation for the time  hr and time step of 0.1 hr. The corresponding seepage values were evaluated as data base of genetic programming. The modelling parameters and their corresponding values are shown in Table 1. Due to the symmetry of the problem and to save time, half of a cavern was modelled. Cavern dimension and initial and boundary conditions for numerical modelling of the gas and oil seepage are shown in Figure 3. Hydrocarbon flow is driven by the difference of Dirichlet pressure boundary conditions which is hydrostatic () for the oil seepage modelling and constant () for the gas seepage modelling and hydrostatic pressure of groundwater. Dirichlet boundary condition, and , is considered for the groundwater table. boundaries are impermeable and given as Neumann no-flow boundary conditions. Rock mass is initially fully saturated with the water and the water pressure in the domain is hydrostatic. The free quad finite elements mesh was used for the modelling. Grids in the area of seepage passage were refined to smaller elements to have more accuracy. Mesh dependency tests were carried out for each case and the meshes eventually used were justified by the quality of the results. Three and four levels of the groundwater level were considered for each of the gas and oil seepage modelling, respectively, where minimum groundwater level was considered to be 2 m above the cavern crown and maximum level is 1 m below the level that no seepage will occur. In order to overestimate the hydrocarbon seepage to have higher factor of safety and for the sake of simplicity, several assumptions were taken into account in the oil and gas seepage modelling as follows:(i)The equivalent-continuum approach modelling was used and no distinction was made between fractures and the matrix blocks and fluids were assumed to flow through the whole system.(ii)The porous medium, representing the rock, was considered homogeneous and isotropic for the oil seepage and homogenous and anisotropic () for the gas seepage modelling.(iii)The rock and the fluids were considered incompressible.(iv)The relative permeabilities and the capillary pressure function of the wetting and nonwetting phases were considered based on Brooks and Corey’s coloration and .(v)The capillary pressure was consumed equal to the entry capillary pressure in the gas seepage modelling and Klinkenberg effect was neglected.(vi)The gases were considered ideal and solubility of the gas in the water and pressure drop of the gas were neglected.(vii)Water density and viscosity were considered 1000 Kg/m3 and 0.001 Pa·m, respectively.The allowable seepage value per unit length (1 m) of the half of the cavern with the stated dimension (Figure 3) would be 0.042 m3/24 h. Therefore, the permeability values were chosen in a domain in which seepage values are close to the allowable seepage value. Corresponding porosity for each permeability value was considered based on Archie’s formulas as follows [19]:where is in millidarcy, mD.

Table 1: Values and domain of the selected parameters.
Figure 3: Scheme of the caverns dimensions and initial and boundary conditions for the numerical modelling of oil (a) and gas (b) seepage.

Since hydrogeologic characteristics of the limestone rocks of Iran are poorly referenced, the entry capillary pressure was measured by function of capillary pressure data in the Edwards formation, Jourdanton field for limestone which is close to the limestone in Iran petrophysically and mineralogically [20]. Each capillary entry pressure value was obtained by the function using a specific permeability and its corresponding value of porosity () as follows [21]:The values of interfacial tension () for the oil-water system and the gas-water system were considered 48 and 50 dyn/cm, respectively [22]. The corresponding values of irreducible water saturation for each porosity value were obtained by Holmes [23] equation: where the maximum value was considered as 0.3. and constant were considered 1 and 0.005–0.06, respectively [23]. Oil and vapor gas density and viscosity and ratio for the limestone were considered based on [2426]. Figures 4 and 5 show the examples of the oil and gas seepage modelling for the parameters presented in Table 2 to give a view of the numerical modelling.

Table 2: Corresponding parameters of the oil and gas seepage modelling example.
Figure 4: Finite element mesh as well as pressure and oil saturation contours of the example of oil seepage modelling.
Figure 5: Pressure and gas saturation contours of the example of gas seepage modelling.

GPTIPS, an open-source MATLAB toolbox for genetic programming (GP) technique, was used to generate prediction equations based on the hydrocarbons seepage value after 24 hr corresponding to each numerical modelling. To reach the stated goal, the whole data set was separated into two, training and testing set (80% and 20% of the whole data set, resp.). The training set was used to evaluate the final or optimum computer programs (CPs) while the testing set was employed to validate the reliability of the GP model. The optimum combination of the values for the set of parameters such as population size, number of generations, function set, mutation rate, crossover rate is achieved by the performance of several trials. The mean absolute percent error (MPE) between the seepage values evaluated by numerical modelling (COMSOL) and the values returned by the GP generated equation was used in the evaluation stage as the fitness function. The mathematical phrase of the best and simplest generated computer program (CP) by GP for the seepage value of the gas and oil was considered as final equation.

6. Results and Discussion

To have more accurate formulations, two equations were presented for the gas seepage where the ground water level is low and the gas seeps from all parts of the gas filled space of the cavern and where the gas seepage is not from all parts of the gas filled section due to high water level above the cavern. Two equations were presented for the oil seepage for two permeability value intervals as well.

Table 3 represents the simplest and the best mathematical phrases generated by the GP where is the intrinsic permeability, in mD, is the pressure difference of oil at its free level and groundwater hydrostatic pressure, in meters of water, is the water level above the cavern crown, in m, is the gas pressure, in bar, is oil density, in g/cm3, is dynamic viscosity of oil, in Pa.m, and is the vertical to horizontal permeability ratio.

Table 3: GP equations of the oil and gas seepage value (m3/24 hr).

Predicted seepage values by the GP equations versus actual values of the seepage for the training and test sets (80 and 20% of data set) as well as relative percent errors of equations   ()–() in Table 3 for the whole data set are shown in Figures 6 and 7, respectively.

Figure 6: Predicted values of the oil and gas seepage (m3/24 h) by the GP equations of ()–() in Table 3 ((a)–(d), resp.) versus actual values for the training and test sets.
Figure 7: Relative percent errors of the equations ()–() in Table 3 ((a)–(d), resp.) for the whole data set.

By considering some parameters (e.g., oil density, viscosity, and gas pressure) as given values and solving the proposed equations for the seepage value less than the allowable one, prejudgment can be done and required permeability value or groundwater level above the cavern can be determined.

The results showed that, in order to control the oil seepage by permeability control technique such as grouting, the permeability of the rock must be less than 100 mD (10−13 m2) with proper water level above the cavern; otherwise the seepage of the stored products would be much more than allowable seepage value. Since grouting cost to have the safe value of permeability is much expensive and complicated, it is better to use hydrodynamic control technique and to locate the cavern deep enough below the groundwater level with a good margin of safety so that no seepage will occur. Since groundwater level decreases due to water leaking to the cavern, it has to be maintained at its original level. Therefore, the cavern must be equipped with the artificial system for supplying water which can be done by injecting water through water curtain systems above the cavern or vertical wells.

7. Conclusion

In this paper prediction equations for the hydrocarbons (oil and gas) seepage from the Iranian unlined rock caverns (URCs) are presented using genetic programming on the basis of the data gathered by numerical modelling of hydrocarbon seepage for the time  hr and different physicochemical properties of the hydrocarbons and hydrogeological properties of the host rock. The seepage control technique can be selected and prejudgment can be done by the application of the presented equations. Since the dimensions of the cavern for numerical modelling (Iranian caverns) are common, the results can be used for crude oil storage unlined rock caverns (URCs) worldwide. The coefficient of determination () and relative percent errors of the proposed equations show their ability in the hydrocarbon seepage predication from URCs in the carbonated rocks with the intrinsic permeability between 10−13 and 10−15 m2 (1–100 mD).


Coefficient of Multiple Determination ()

The coefficient of multiple determination () measures how successful the fit is in explaining the variation of data. Thus the closer is to unity indicates the better model performance and fit.where is the number of data, and are the actual and predicted values, respectively, and is the average of ().

The adjusted , the degree of freedom, is a modified version of that has been adjusted for the number of predictors in the model and is generally the best indicator of the fit quality when you add additional coefficients to a model. It is always lower than or equal to .where is the number of predictors.

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.


The authors would like to acknowledge the financial support from Iranian Oil Terminals Company and thank Professor Seyed Majid Hassanizaeh from Utrecht University for thoughtful comments and helpful discussions.


  1. C.-O. Morfeldt, “Storage of petroleum products in man-made caverns in Sweden,” Bulletin of the International Association of Engineering Geology, vol. 28, no. 1, pp. 17–30, 1983. View at Publisher · View at Google Scholar · View at Scopus
  2. S. M. Haug, “Storage of oil and gas in rock caverns below the groundwater table- general design development,” in Underground Constructions for the Norwegian Oil and Gas Industry, pp. 19–25, Norwegian Tunnelling Assossiation, 2007. View at Google Scholar
  3. B. Aberg, “Prevention of gas leakage from unlined reservoirs in rock,” in Rockstore, vol. 77, pp. 399–413, Pergamon Press, 1977. View at Google Scholar
  4. N. O. Midtlien, “Cavern storage excavtion-sture,” in Underground Constructions for the Norwegian Oil and Gas Industry, pp. 47–56, Norwegian Tunnelling Assossiation, 2007. View at Google Scholar
  5. H. Liu, J. C. Lee, and B. R. Li, “High precision pressure control of a pneumatic chamber using a hybrid fuzzy PID controller,” International Journal of Precision Engineering and Manufacturing, vol. 7, pp. 7–13, 2007. View at Google Scholar
  6. C. Yu, S. C. Deng, H. B. Li, J. C. Li, and X. Xia, “The anisotropic seepage analysis of water-sealed underground oil storage caverns,” Tunnelling and Underground Space Technology, vol. 38, pp. 26–37, 2013. View at Publisher · View at Google Scholar · View at Scopus
  7. C. Y. Wang and C. Beckermann, “A two-phase mixture model of liquid-gas flow and heat transfer in capillary porous media-I. Formulation,” International Journal of Heat and Mass Transfer, vol. 36, no. 11, pp. 2747–2758, 1993. View at Publisher · View at Google Scholar · View at Scopus
  8. M. K. Zahoor, M. N. Derahman, and M. H. Yunan, “Implementation of Brooks and Corey correlation in water wet case with immobile wetting phase,” NAFTA, vol. 60, pp. 435–437, 2009. View at Google Scholar
  9. K. Li, Theoretical Development of the Brooks-Corey Capillary Pressure Model from Fractal Modeling of Porous Media, SPE, 2004.
  10. O. B. Baker, H. W. Yarranton, and J. Jensen, Practical Reservoir Engineering and Characterization, Gulf Professional Publishing, 1st edition, 2015.
  11. J. Rutqvist, Y.-S. Wu, C.-F. Tsang, and G. Bodvarsson, “A modeling approach for analysis of coupled multiphase fluid flow, heat transfer, and deformation in fractured porous rock,” International Journal of Rock Mechanics and Mining Sciences, vol. 39, no. 4, pp. 429–442, 2002. View at Publisher · View at Google Scholar · View at Scopus
  12. V. Reichenberger, H. Jakobs, P. Bastian, and R. Helmig, “A mixed-dimensional finite volume method for two-phase flow in fractured porous media,” Advances in Water Resources, vol. 29, no. 7, pp. 1020–1036, 2006. View at Publisher · View at Google Scholar · View at Scopus
  13. M. Diaz-Viera, “COMSOL implementation of a multiphase fluid flow model in porous media,” in Proceeding of the COMSOL Conference, 2008.
  14. M. F. El-Amin and S. Shuyu, “Effects of gravity and inlet/outlet location on a two-phase cocurrent imbibition in porous media,” Journal of Applied Mathematics, vol. 2011, Article ID 673523, 2011. View at Publisher · View at Google Scholar · View at MathSciNet
  15. D. Droste, F. Lindner, C. Mundt, and M. Pfitzner, “Numerical computation of two-phase flow in porous media,” in Proceeding of the COMSOL Conference, pp. 1–6, 2013.
  16. Y. Cao, E. Birgitte, and R. Helmig, “Fractional flow formulation for two phase flow in porous media,” in International Research Training Group GRK 1398/1, Non-Linearities and Upscaling in Porous Media, 2007.
  17. J. R. Koza, Genetic Programming: on the Programming of Computers by Means of Natural Selection, MIT Press, Cambridge, Mass, USA, 1976. View at Publisher · View at Google Scholar · View at Scopus
  18. E. Ghotbi Ravandi, R. Rahmannejad, A. E. Feili Monfared, and E. Ghotbi Ravandi, “Application of numerical modeling and genetic programming to estimate rock mass modulus of deformation,” International Journal of Mining Science and Technology, vol. 23, no. 5, pp. 733–737, 2013. View at Publisher · View at Google Scholar · View at Scopus
  19. S. Gholinezhad and M. Masihi, “A physical-based model of permeability/porosity relationship for the rock data of Iran southern carbonate reservoirs,” Iranian Journal of Oil & Gas Science and Technology, vol. 1, pp. 25–36, 2012. View at Google Scholar
  20. J. W. Amyx, D. M. Bass, and R. L. Whiting, Petroleum Reservoir Engineering Physical Properties, McGraw-Hill, New York, NY, USA, 1960.
  21. T. Ahmed, Reservoir Engineering Handbook, Elsevier, 4th edition, 2010.
  22. J. H. Schon, “Physical Properties of Rocks,” in Handbook of Petroleum Exploration and Production, Elsevier B. V., Amsterdam, The Netherlands, 2011. View at Google Scholar
  23. M. Holmes, A. Holmes, and D. Holmes, “Relationship between porosity and water saturation: methodology to distinguish mobile from capillary bound water two different rock types,” American Association of Petroleum Geologists (AAPG), pp. 1–11, 2009. View at Google Scholar
  24. M. Sattarin, H. Modarresi, and M. Teymori, “New viscosity correlations for dead crude oils,” Petroleum & Coal, vol. 49, pp. 33–39, 2007. View at Google Scholar
  25. N. P. Cheremisinoff and A. Davletshin, A Guide to Safe Material and Chemical Handling, Jon Wiley & Sons, 2010.
  26. P. A. Domenico and F. W. Schwartz, Physical and Chemical Hydrogeology, John Wiley & Sons, 1990.