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

Authors: Ebrahim Ghotbi Ravandi; Reza Rahmannejad; Saeed Karimi-Nasab; Amir Sarrafi
Journal: Geofluids
 (2017)

---

## Abstract

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 (R2) 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.


---

## Body

## 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 https://doi.org/10.1155/2017/6803294). 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 u are denoted without index. With introduction of the mixture quantities (see [7]) the conservation of mass with the porosity (ϕ) is defined as(1)∂ϕρ∂t+∇·ρu=0.Conservation of momentum using Darcy velocity is as follows in which the dynamic viscosity μ and the pressure p are also mixture quantities (see [7]):(2)u=-Kμ∇p-ρkg∇z1μ=swkrwμw+snwkrnwμnwρk=ρwkrwμw+ρnwkrnwμnw∇p=krwμw∇pw+krwμw∇pnw,where ρk is the kinetic mixture density, K is the intrinsic permeability, s is the saturation and the subscripts of w and nw are related to the wetting and nonwetting phase, respectively, kr is the relative permeability, g is the gravitational acceleration vector, and z is the depth.The diffusive mass fluxj connects the mixture mean velocity with the velocity of the individual phases:(3)ρwuw=krwμwρu+jρnwunw=krnwμnwρu-j.Wang and Beckerman [7] introduced the diffusive flux (j) as follows:(4)j=DC∇sw+fswKΔρνnwg,where ν is the kinematic viscosity, f is the hindrance function for phase migration and separation, and DC is the capillary diffusion coefficient as a function of the wetting phase saturation:(5)DC=krwμwKsw-1∂pc∂sw,where pc is capillary pressure. The transport equation is as follows, where cw is the fluid content of the wetting phase:(6)∂ϕcw∂t+∇·cwu=∇·DC∇cwcw=ρwsw.Several scientists have tried to derive a functional correlation for the relative permeability and the capillary pressure (pc) 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 (pec) and the wetting phase saturation that empirically has been found to be appropriate for the drainage process as follows [8, 9]:(7)pc=pecsw-1/λkrw=sw3+2/λkrnw=1-sw21-sw2+λ/λ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 [11–15]. 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 Figure1. 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 T=2000 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 (Sw) fronts and contour lines for the time T=2000 days and the time step (Δt) 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 T=24 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 (Pnw=Pgas+ρoilgh) for the oil seepage modelling and constant (Pgas) for the gas seepage modelling and hydrostatic pressure of groundwater. Dirichlet boundary condition, Pw=0 and Sw=1, is considered for the groundwater table. Γnoflow 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 (Kx=Ky≥Kz) 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λ=2.(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]:(8)K=2.5510ϕ5.65K=9.3510ϕ5.65,where K is in millidarcy, mD.Table 1
Values and domain of the selected parameters.
Intrinsic permeability (m2)
Porosity
Irreducible water saturation
Nonwet phase density (Kg/m3)
Nonwet phase viscosity (Pa⋅s)
Gas pressure (bar)
K
z
/
K
x
Oil
1
e
-
13
0.15–0.2
0.05–0.3(3)∗
800–950(3)
0.1–5.42
—
—
5
e
-
14
0.12–0.17
0.05–0.3(3)
800–950(3)
0.1–5.42
—
—
1
e
-
14
0.1–0.15
0.05–0.3(3)
800–950(3)
0.1–5.42
—
—
5
e
-
15
0.08–012
0.05–0.3(3)
800–950(3)
0.1–5.42
—
—
1
e
-
15
0.05–0.1
0.1–0.3(3)
800–950(3)
0.1–5.42
—
Gas
1
e
-
14
0.1–0.15
0.05–0.3
3.125
2
e
-
5
1-2(3)
0.5–1(3)
5
e
-
15
0.08–0.12
0.05–0.3
3.125
2
e
-
5
1-2(3)
0.5–1(3)
1
e
-
15
0.05–0.1
0.1–0.3
3.125
2
e
-
5
1-2(3)
0.5–1(3)
∗The number in the parentheses represents the number of levels for the parameter.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 byJ 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 J function using a specific permeability and its corresponding value of porosity (ϕ) as follows [21]:(9)JSw=0.21645PcσKϕ.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: (10)ϕQ×Swc=constant,where the maximum value was considered as 0.3. Q and constant were considered 1 and 0.005–0.06, respectively [23]. Oil and vapor gas density and viscosity and Kz/Kx ratio for the limestone were considered based on [24–26]. 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.
Parameter
Value
Gas pressure (Pgas)
2 bar
Water level above the oil level
6 m
Porosity
0.15
Residual water saturation
0.05
Oil density
950 Kg/m3
Oil viscosity
5.42 Pa⋅m
Gas density
3.125 Kg/m3
Gas viscosity
2
e
-
5 Pa⋅mFigure 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.Table3 represents the simplest and the best mathematical phrases generated by the GP where K is the intrinsic permeability, in mD, Pd is the pressure difference of oil at its free level and groundwater hydrostatic pressure, in meters of water, wl is the water level above the cavern crown, in m, Pgas is the gas pressure, in bar, ρ is oil density, in g/cm3, μ is dynamic viscosity of oil, in Pa.m, and Kz/Kx is the vertical to horizontal permeability ratio.Table 3
GP equations of the oil and gas seepage value (m3/24 hr).
Equation number
Description
Mathematical phrase
Mean percent error
Max percent error
Coefficient of determination (R2)∗
AdjustedR2 (adj_R2)∗
(
∗
)
Oil seepageK = 100–10 md
2.608
K
3
1
0
8
+
-
3.847
P
d
2
1
0
8
-
7.176
1
0
6
K
2
+
0.0001643
P
d
-
7.176
S
w
c
1
0
6
+
0.0001391
log
⁡
K
ϕ
ρ
oil
+
0.0004071
S
w
c
P
d
sin
⁡
ρ
oil
+
2.608
1
0
8
K
+
0.0005122
P
d
-
0.001781
3.32
17
99.4
99.2
(
∗
∗
)
Oil seepageK = 10–1 mD
0.0003502
S
w
c
K
P
d
+
1.961
+
6.479
K
K
P
d
+
1.961
1
0
6
-
3.49
K
+
8.366
/
1
0
5
sin
⁡
sin
⁡
K
-
4.428
exp
⁡
ρ
oil
K
-
K
P
d
10
5
-
2.206
K
+
16.55
/
1
0
6
sin
⁡
μ
oil
+
0.0006279
2.2
13.2
99.8
99.7
(
∗
∗
∗
)
Gas seepage,from all parts
0.09853
w
l
-
0.0878
P
gas
K
-
0.003778
arcsin
⁡
K
z
/
K
x
log
⁡
w
l
+
0.0005139
P
gas
S
w
c
log
⁡
P
gas
K
+
0.002057
-
3.782
ϕ
exp
⁡
P
gas
-
0.0669
P
gas
w
l
2
5.9
18
98
96
(
∗
∗
∗
∗
)
Gas seepage,not from all parts
0.0001787
K
cosh
⁡
w
l
+
0.001714
K
P
gas
-
w
l
cos
⁡
ϕ
-
0.1371
ϕ
tanh
⁡
tanh
⁡
K
z
/
K
x
P
gas
+
0.0002673
K
w
l
K
z
/
K
x
+
0.0076
K
P
gas
-
5.055
1
0
5
12
29
92
90
R
2
∗ and adj_R2 are described in Appendix.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.
(a)
(b)
(c)
(d)Figure 7
Relative percent errors of the equations (∗)–(∗∗∗∗) in Table 3 ((a)–(d), resp.) for the whole data set.
(a)
(b)
(c)
(d)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 timeT=24 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 (R2) 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).



---

*Source: 6803294-2017-06-14.xml*