Document Type: Research Paper
Authors
^{1} Assistant Professor, Faculty of Chemical Engineering, Tarbiat Modares University, Tehran, Iran.
^{2} MS Graduate Student, Faculty of Chemical Engineering, Tarbiat Modares University, Tehran, Iran.
^{3} PhD Candidate, Faculty of Chemical Engineering, Tarbiat Modares University, Tehran, Iran.
^{4} Natural Gas Storage Company, Tehran, Iran.
Abstract
Keywords
Main Subjects
The general process of natural gas storage in underground reservoirs involves injecting natural gas in reservoirs during the periods of reduced demand and extracting it in rising demand periods (Tek, 1989; Dietert & Pursell, 2000; Aminu et al., 2017). Demand fluctuations arising from climatic conditions have raised the natural gas storage in underground reservoirs as an economically efficient method of gas storage (Gluyas & Mathias, 2014; Van der Meer, 2005). Generally, the methods of natural gas storage include storing natural gas in depleted oil and gas reservoirs, storage in salt domes, and storage in aquifers (Ji et al., 2017; Lorenz & Müller, 2003; Michael et al., 2010; Shahmorad et al., 2016). Among these methods, natural gas storage in aquifers is very important due to the suitable storage conditions and the ability of stored natural gas trapping (Kvamme et al., 2007; Reidel et al., 2002; Kumar et al., 2005). Aquifers because of the high volume for storage and the presence of water to provide pressure and prevent the leakage of injected gas are one of the most suitable options for natural gas storage in underground reservoirs (Thompson et al., 2009; Bruant et al., 2002). It should be mentioned that natural gas storage in an aquifer core sample has not been studied yet.
The effects of some geological parameters such as permeability, formation dip angle, porosity, formation thickness, and other parameters were investigated in previous studies (Wang & Holditch, 2005), and the qualification of the several reservoirs was checked as a place for underground gas storage processes (Muonagor & Anyadiegwu, 2014; Coffin & Lebas, 2008; Khashayar, 2010). The results show that all parameters can affect the process, but the permeability and thickness of formation are dominant parameters influencing gas injection (Wang & Holditch, 2005). In other words, recovered gas volume percentage and cumulative gas injection in gas storage process are strongly linked to primary reservoir parameters (Wang & Holditch, 2005). In addition, partially depleted reservoirs, due to their high field pressure, are a good candidate to store gas (Soroush & Alizadeh, 2007; Azin et al., 2008).
However, conducted studies on the natural gas storage in aquifers are limited, and many aspects of this process have not been examined yet. Hence, the aims of this work are the simulation of natural gas storage in an aquifer through core scale simulation by Eclipse 300 software and the investigation of the parameters affecting this process through the design of experiments (DOE) for the first time.
DOE is one of the best statistical and mathematical methodologies helping researchers to plan minimum required tests for evaluating and screening the effective parameters and their interactions on a specific process (Gharibshahi et al., 2015; Antony, 2014). Also, this method can be used to optimize the process and find the best combination of effective factors to reach a suitable system response (Montgomery, 2017).
Reducing the time and cost, studying the complex and difficult conditions of a process, and obtaining an accurate and complete form of the system are the advantages of simulations compared to experimental works (Jafari, 2008; Gharibshahi et al., 2015). Therefore, in this study, Eclipse 300, useful software of the petroleum industry, was utilized to simulate natural gas storage in a core sample considered as an aquifer. Moreover, for validating the accuracy of simulation, the numerical results were compared with experimental data.
Additionally, in this research, all of the parameters influencing the process were examined and screened by using PlackettBurman method. Finally, two models were developed using response surface method (RSM) for the statistical analysis and optimization of the most significant factors affecting the stored gas volume and recovered gas volume as the responses of the system.
2. Numerical implementation
2.1. Geometry
In order to create the geometry, a horizontal core plug was simulated using Eclipse 300. The porous media is a core with a diameter of 38 mm and a length of 150 mm. Furthermore, the permeability of the core is 172 md, and its porosity is about 20%. This geometry was the same as the porous media that Junho Oh et al. (2013) used in their experimental work. Figure 1 shows the schematic of the porous media used in this paper.
Figure 1
A schematic of the porous medium.
2.2. Grid independency check
The number of elements in gridding the model can significantly affect the numerical results (Crochet et al., 2012). For this reason, the number of elements must be enough to be able to show the changes during the simulation well. By increasing the number of elements in the model, the accuracy of the numerical results increases, so the time of calculation rises too (Mousavi et al., 2009). Hence, checking the grid independency of the model is a necessary work in every simulation process. Therefore, in order to investigate the sensitivity of numerical results in comparison with the number of elements, the number of cells in xdirection is changed, and three models with different number of cells were generated. Then, the numerical results of CO_{2 }injection into the porous media were determined at a constant flow rate (5 cc/min) and in the condition of closed production wells. Table 1 shows the results of gas saturation in the porous media as the main parameter to check the grid independency.
Table 1
Results of the grid independency test.
Number of elements in xdirection 
Number of elements in ydirection 
Number of elements in zdirection 
Number of Cells 
Average pressure (atm) 
Gas saturation (%) 
40 
15 
15 
9000 
302.5544 
8.1147 
50 
15 
15 
11250 
305.8419 
8.1712 
60 
15 
15 
13500 
305.8529 
8.1727 
According to Table 1, it can be seen that the variations of gas saturation and the core pressure after 11250 elements are very small. Therefore, to save the calculation time and computational cost, grid 2 with 11250 cells (50 elements in xdirection and 15 elements in y and zdirections) was chosen as the main grid. Figure 2 shows the discretized model produced in Eclipse 300 software.
To have a condition same as the experimental work, all the surface area of the left side of the core was assumed as the inlet port. Because of the existence of 157 elements in the left side of the core, all of these elements were assumed as an injection port. Therefore, the injection flow rate of the inlet port is the sum of the all flow rates in each element. In addition, all the surface area of the right side of the core was assumed as the outlet port at which the fluid is produced at a constant pressure of 98.96 atm.
Figure 2
A schematic of the discretized model.
2.3. The governing equations of twophase flow of water and gas
The general governing equations of multiphase flow into the system (gas flow in a saturated water porous medium) based on the conservation of mass law were calculated as follows (Soave, 1972):
(1) 

(2) 
In these equations, the subscript w shows the properties of water, and subscript g shows the properties of gas. Also, ρ is density, and μ is viscosity; K, P, P_{c}, s, and q stand for permeability, pressure, capillary pressure, saturation, and flow rate respectively. Moreover φ, z, g, and t represent porosity, height, gravity, and time respectively.
To solve these equations, according to the discretized 3D porous media, the full implicit (FULL IMPES) method was used to simultaneously solve these equations with the equation of state (EOS). In these equations, the pressure of gas and water were related to each other by capillary pressure (P_{c}) as given in Equation 3:
(3) 
2.4. Equation of state (EOS)
Cubic equations of state were used to determine the thermodynamic properties and phase equilibrium calculations for compounds containing heavy and complex hydrocarbons widely (Ahmed, 2006; Soave, 1972). Thus, in this research, PengRobinson EOS is used to determine the thermodynamic properties of the existing fluids (natural gas and water) in the system.
2.5. Available fluids in the system
In this work, for the validation of the numerical results with experimental data, CO_{2} was assumed as the injected fluid, and for the other tests, i.e. screening and designing the experiment tests, natural gas was assumed as the injected fluid. Similar to the experimental work, the simulation was conducted at a temperature of 40 °C and at ambient pressure, and the properties of available fluids in this study were determined in the above conditions. The physical properties of the available fluids in the system are shown in Tables 2 and 3.
Table 2
Composition of natural gas.
Component 
Mol. % 
Methane 
90.490 
Ethane 
6.707 
Propane 
0.359 
Nitrogen 
1.027 
Carbon dioxide 
1.417 
Table 3
Physical properties of the available fluids in the system at 40 °C.
Value 
Property 
Nanoparticle 
95.21 
Density, ρ (kg/m^{3}) 
Natural gas 
15.68×10^{6} 
Viscosity, μ (Pa.s) 

17.6 
Molecular weight, Mw (kg/kmol) 

1.697 
Density, ρ (kg/m^{3}) 
CO_{2} 
15.65×10^{6} 
Viscosity, μ (Pa.s) 

44.01 
Molecular weight, Mw (kg/kmol) 

992.3 
Density, ρ (kg/m^{3}) 
Water 
0.000653 
Viscosity, μ (Pa.s) 

18.05 
Molecular weight, Mw (kg/kmol) 
2.6. Validation of numerical results
To validate the accuracy of the numerical results, after checking the grid independency, CO_{2 }injection into theporous medium was carried out in the conditions similar to the experimental work. To this end, the flow rate of injection was considered equal to 5 cc/min. Due to the presence of 157 injection wells, the injection flow rate of each well was calculated to be 0.0318 cc/min. After the injection of 4 pore volumes (PV) of CO_{2 }into the porous medium, the results of gas saturation into the porous medium were compared with the experimental data as shown in Figure 3. It can be seen from Figure 3 that the numerical results are in good agreement with the experimental data, and the maximum relative error between them is about 9.5 %. Also, the relative error between the numerical results and the experimental data is about 5% throughout the core. Figure 4 shows the volume fraction contour of gas saturation in the porous media. It can be seen that at the beginning of the core, water was well displaced by the gas. Despite drastic changes of saturation in the length of the core, because of the small size of the core in the vertical direction, significant changes of saturation is not observed in this direction. It is worth mentioning that the residual water saturation in core is about 40%, which shows that a lot of water was displaced by the injected gas at the beginning of the core.
Figure 3
Comparison of the numerical results with the experimental data.
Figure 4
The volume fraction contour of gas saturation in the porous medium.
3. Design of experiments and investigation of the effective parameters
Design of experiments is a statistical and mathematical method that can be used to plan minimum requirement tests to investigate the effect of some parameters and their interactions in a process and to optimize them to reach the best condition in a system (Gharibshahi et al., 2015). Therefore, in order to investigate the effect of various parameters on the natural gas storage process, two DOE methods were used herein. At first, some parameters affecting the process were determined. By considering different criteria presented in the literature related to the natural gas storage processes, seven parameters, including the injected flow rate, irreducible water saturation (S_{Wirr}), pressure, salinity, temperature, porosity, and permeability were determined. All of these parameters and their levels are tabulated in Table 4. It should be mentioned that the range of the pressure, salinity, and temperature are selected according to one of Iranian aquifer, and the range of the other parameters are chosen regarding the simulated core properties and operational feasibility. Because of the large number of effective parameters in this process, which increases the number of requirement tests in the DOE plan, PlackettBurman method was first used for screening parameters and finding the most important parameters in this process to make a good DOE plan and study this process well.
Table 4
Effective parameters and their levels.
No. 
Parameter 
Level 1 
Level 2 
Level 3 
Level 4 
Level 5 
1 
Flow rate (cc/min) 
0.1 
0.3 
0.5 
0.7 
0.9 
2 
S_{wirr} 
0.1 
0.2 
0.3 
0.4 
0.5 
3 
Pressure (psi) 
1000 
1500 
2000 
2500 
3000 
4 
Salinity (ppm) 
100000 
125000 
150000 
175000 
200000 
5 
Temperature (°C) 
25 
36 
47 
58 
69 
6 
Porosity 
4 
8 
12 
16 
20 
7 
Permeability (md) 
0.2 
0.6 
1 
1.4 
1.8 
4. Results and discussion
4.1. Parameter screening
The PlackettBurman method is one of the DOE methods widely used for screening the effective parameter in a process (Montgomery, 2017). Due to the presence of seven parameters in this process, the most effective parameters can be found by eight tests suggested by this method, and we will have a good DOE plan with a fewer number of tests to save time and cost. In this method, two levels must be considered for each parameter (the upper and lower levels of each parameter). Table 5 shows the tests suggested by this method and the level of each factor. In this table, the upper level of each parameter is shown by +1, and the lower level of each parameter is represented by 1.
In order to investigate the effect of parameters, an objective function must be considered as the main criteria for the comparison of the parameters and their interactions. In this study, storage gas volume (gas saturation) in core is considered as the response of the system. This saturation is determined after injecting natural gas at a constant flow rate of 5 cc/min into the core for 30 min. The results of gas saturation in each test are shown in Table 6.
Table 5
PlackettBurman design developed for screening the effective parameters.
Run 
Flow rate (cc/min) 
S_{wirr} 
Pressure (psi) 
Salinity (ppm) 
Temperature (°C) 
Porosity 
Permeability (md) 
1 
+1 
+1 
+1 
+1 
+1 
+1 
+1 
2 
+1 
+1 
1 
1 
1 
1 
+1 
3 
+1 
1 
+1 
+1 
1 
1 
1 
4 
+1 
1 
1 
1 
+1 
+1 
1 
5 
1 
+1 
+1 
1 
+1 
1 
1 
6 
1 
+1 
1 
+1 
1 
+1 
1 
7 
1 
1 
+1 
1 
1 
+1 
+1 
8 
1 
1 
1 
+1 
+1 
1 
+1 
Table 6
Results of gas saturation in the PlackettBurman tests.
Test 
1 
2 
3 
4 
5 
6 
7 
8 
Gas saturation 
0.1461 
0.1368 
0.2114 
0.2100 
0.0513 
0.0561 
0.1152 
0.1201 
After performing the tests and recording the results, the effect of each parameter can be found by Equation 4:
(4) 
The amount of effect for each parameter 
where, n is the number of tests, and (y+) and (y) are the response of the system when the desired parameter is at the maximum and minimum level respectively. The absolute value of value obtained for each parameter indicates the importance of each parameter. Table 7 lists the amount of effect of each parameter on the response of system in the PlackettBurman method.
Table 7
The amount of effect for each parameter in PlackettBurman method.
Parameter 
Flow rate (cc/min) 
S_{wirr} 
Pressure (psi) 
Salinity (ppm) 
Temperature (°C) 
Porosity 
Permeability (md) 
Amount of Effect 
0.0902 
0.0667 
0.00902 
0.00091 
0.00027 
0.0016 
0.0035 
As listed in Table 7, by considering the absolute value of the obtained values, the flow rate of the injected fluid has the greatest impact on the amount of stored gas. The first reason for this is that by increasing the flow rate, the value of capillary number increases, and water is displaced by the gas more efficiently and exits from the core. The second reason in this case is injecting more volume of gas during a certain time at higher flow rates. By increasing the volume of the injected gas in the core, a larger amount of water will be displaced by the gas. Also, increasing the irreducible water saturation makes a reduction in gas saturation in the core. This phenomenon occurs also due to increasing the capillary pressure exponentially when water saturation is close to its final quantity. By increasing the capillary pressure, water displacement becomes more difficult.
Furthermore, the amount of stored gas in the core decreases by increasing the pressure. This may happen due to the increase in gas density by increasing the pressure and decreasing its mobility. Increasing the salinity of the water reduces the amount of stored gas too because by increasing the salinity of water, the viscosity of water increases, which makes a reduction in water mobility. Temperature has a positive effect on the gas saturation, and, by increasing the temperature, gas saturation rises, which is due to the improvement of water viscosity by raising the temperature. In the same conditions of gas injection, the higher the pore volume of the core is, the higher the volume of the existing water and water saturation becomes. As a result, by increasing the porosity of the core, gas saturation drops. Moreover, since by increasing the permeability of the core, the mobility of both phases rises, the permeability of the core has a positive impact on gas saturation.
Finally, according to Table 7, it can be found out that the flow rate of injection fluid and the irreducible water saturation have the greatest impact on the process, and porosity, temperature, and salinity can be ignored. It is worth mentioning that although pressure and permeability have little effect on the process, these two parameters are considered in the final DOE plan using the response surface method (RSM).
4.2. Parametric study and optimization results
After screening the parameters affecting the process and finding the most significant parameters, four parameters, including flow rate, irreducible water saturation, pressure, and permeability were considered to plan a systematic study on natural gas storage. To this end, response surface method was selected for the design of the experiments. The purpose of using this method is the optimization of the process and finding the optimum and the best level of each parameter to reach the maximum gas storage.
RSM is a statistical method that investigates the relationship between several parameters and one or more responses introduced by G. E. P. Box and K. B. Wilson in 1951. In this method, the range of numerical parameters is considered as a continuous interval, and the effect of all integers in the interval is examined in the DOE plan (Yu & Sepehrnoori, 2013; Eghbali et al., 2013).
According to the selected parameters, to design the required simulation tests, Design Expert 7 software was used. Table 8 shows the tests designed by the RSM method. In this study, two parameters, including stored gas volume (the volume (cc) of gas stored in the porous core) and the recovered gas volume (the ratio (%) of the produced gas volume to the volume of the stored gas in the core were considered as the responses of the system. According to the suggested conditions for each test in Table 8, all of the tests were carried out by Eclipse 300 software, and the relevant stored gas volume and the recovered gas volume after 670 s of gas injection were determined as the corresponding responses. All of the relevant responses for each test are tabulated in Table 8.
Table 8
RSM designed plan developed for the systematic study of natural gas storage.
No. 
Flow rate (cc/min) 
S_{wirr} 
Pressure (psi) 
Permeability (md) 
Stored gas volume (cc) 
Recovered gas volume (%) 
1 
0.3 
0.4 
2500 
6 
5.11277 
42.0102 
2 
0.3 
0.2 
2500 
0.6 
5.24595 
40.5805 
3 
0.1 
0.3 
2000 
1 
2.16283 
26.5618 
4 
0.5 
0.3 
2000 
1.8 
8.51515 
43.4831 
5 
0.3 
0.2 
1500 
1.4 
5.60272 
37.5109 
6 
0.7 
0.2 
1500 
1.4 
11.221 
47.2638 
7 
0.7 
0.2 
2500 
0.6 
10.7287 
48.533 
8 
0.5 
0.1 
2000 
1 
8.41131 
45.4043 
9 
0.5 
0.3 
2000 
1 
8.19776 
45.6909 
10 
0.5 
0.3 
2000 
0.2 
7.49102 
50.6018 
11 
0.5 
0.3 
2000 
1 
8.19776 
45.6909 
12 
0.3 
0.2 
2500 
1.4 
5.59593 
36.9012 
13 
0.3 
0.4 
2500 
1.4 
5.46574 
39.1357 
14 
0.5 
0.3 
2000 
1 
5.44574 
40.7779 
15 
0.5 
0.3 
3000 
1 
5.33926 
37.5281 
16 
0.7 
0.4 
2500 
0.6 
10.4053 
50.03555 
17 
0.7 
0.4 
2500 
1.4 
10.7721 
49.7757 
18 
0.5 
0.5 
2000 
1 
7.88055 
53.9742 
19 
0.3 
0.2 
1500 
0.6 
5.29993 
41.2858 
20 
0.5 
0.3 
1000 
1 
8.18468 
45.2505 
21 
0.7 
0.2 
1500 
0.6 
10.8064 
49.3634 
22 
0.7 
0.4 
1500 
0.6 
10.4715 
50.8601 
23 
0.3 
0.4 
1500 
1.4 
5.50615 
39.9188 
24 
0.5 
0.3 
2000 
1 
8.19776 
45.6909 
25 
0.9 
0.3 
2000 
1 
13.5873 
50.931 
26 
0.7 
0.2 
2500 
1.4 
11.1569 
47.1779 
27 
0.7 
0.4 
1500 
1.4 
10.8349 
49.1503 
28 
0.3 
0.4 
1500 
0.6 
5.16277 
42.7484 
4.3. Analysis of variance
The analysis of variance is a method mainly used to determine the important parameters influencing the responses through the analysis of the total variance to the right components and measuring the relative effect of each of them (Mason et al., 2003). The results of the analysis of variance for two responses are shown in Tables 9 and 10. These results are obtained based on the backward method in which insignificant terms are automatically ignored. Also, the quadratic statistical model is considered in the model presented by the software. By choosing the quadratic model, the interactions between two different parameters like AB, BC etc. and similar parameters such as A^{2}, B^{2} etc. are considered in the predicted model.
In these tables, the degree of freedom of each parameter is defined as follows:
(5) 
where, N is the number of levels for each factor. The sum of squares for each factor is calculated by Equation 6:
(6) 
where, is the standard average of the results at level j for each parameter, and represents the total standard average. The variance and the percent contribution of each parameter are also calculated by Equations 7 and 8:
(7) 

(8) 
The Fvalue is the ratio of the mean square error to residual error and is usually used to determine the significance of a parameter. When the Fvalue of a parameter increases, the importance of the parameter rises. This entry is vice versa for the pvalue, so by decreasing the pvalue, the importance of that parameter increases. According to Tables 9 and 10, it can be seen that the flow rate has the greatest impact on the stored gas volume and recovered gas volume, and the other three parameters, namely pressure, permeability, and irreducible water saturation, have the least impact on the process.
Table 9
ANOVA analysis of stored gas volume.
Source 
Sum of squares 
DOF 
Mean square 
Fvalue 
pvalue 
Model 
92.00 
14 
6.57 
19.21 
< 0.0001 
A Flow rate 
89.03 
1 
89.03 
260.25 
< 0.0001 
BIrreducible water saturation 
0.18 
1 
0.18 
0.53 
0.4798 
C Pressure 
0.76 
1 
0.76 
2.22 
0.1605 
D Permeability 
0.50 
1 
0.5 
1.46 
0.2477 
A^{2} 
0.29 
1 
0.29 
1.84 
0.3775 
B^{2} 
0.59 
1 
0.59 
1.71 
0.2130 
C^{2} 
0.17 
1 
0.17 
0.51 
0.4888 
D^{2} 
0.41 
1 
0.41 
1.21 
0.2909 
Residual 
4.45 
13 
0.34 


96.45 
27 



Table 10
ANOVA analysis of recovered gas volume.
Source 
Sum of squares 
DOF 
Mean square 
Fvalue 
pvalue 
Model 
868.36 
14 
62.03 
14.07 
< 0.0001 
A Flow rate 
608.09 
1 
608.09 
137.9 
< 0.0001 
BIrreducible water saturation 
49.09 
1 
43.09 
9.77 
0.0080 
C Pressure 
15.68 
1 
15.68 
3.56 
0.0819 
D Permeability 
44.88 
1 
44.88 
10.18 
0.0071 
AD 
3.74 
1 
3.74 
1.15 
0.3739 
A^{2} 
45.66 
1 
45.66 
14.05 
0.0067 
B^{2} 
44.16 
1 
44.16 
13.59 
0.0075 
C^{2} 
12.39 
1 
12.39 
3.81 
0.1176 
D^{2} 
11.58 
1 
11.58 
3.56 
0.1290 
Residual 
57.32 
13 
4.41 


925.68 
27 



The two models presented by software based on the actual values of the two responses, namely the stored gas volume and the recovered gas volume, are tabulated in Table 11.
Table 11
Statistical results of the ANOVA for the reduced quadratic model.
Response 
Final equation in terms of coded factors 
pvalue 
R^{2} 
Adjusted R^{2} 
Adequate precision 
CV (%) 
Stored gas volume 
5.24+1.93A0.087B0.18C+0.14D+0.11A^{2}+0.16B^{2}0.085C^{2}+0.13D^{2} 
< 0.0001 
0.9539 
0.9042 
17.996 
10.62 
Recovered gas volume 
44.64+5.03A+1.34B0.81C1.37D+0.048AD1.38A^{2}+1.36B^{2}0.72C^{2}+0.69D^{2} 
< 0.0001 
0.9381 
0.8714 
15.413 
4.73 
The Rsquare values of these two models are 0.9539 and 0.9381. It is obvious that Rsquare values closer to 1 means a more accurate model. Thus, the high values of the predicted and adjusted Rsquare in this paper indicate that the reduced quadratic model is a good representation of the system.
Figures 5 to 7 show the effect of irreducible water saturation, permeability, and pressure versus the flow rate of the injected gas on the stored gas volume respectively. These figures indicate that by increasing the flow rate of the injected gas, at any values of the other parameters, the value of the stored gas volume increases. In other words, increasing the amount of the other parameters has a negligible effect on the stored gas volume compared to the flow rate of the injected gas, which validates the results of the analysis of variance.
Figure 5
The effect of flow rate variation with irreducible water saturation on the amount of stored gas volume.
Figure 6
The effect of flow rate variation with permeability on the amount of stored gas volume.
Figure 7
The effect of flow rate variation with pressure on the amount of stored gas volume.
Figures 8 to 10 depict the effect of irreducible water saturation, permeability, and pressure versus the flow rate of the injected gas on the recovered gas volume respectively. As can be seen, similar to the stored gas volume, the flow rate of the injected gas has a significant impact on the recovered gas volume, and the other parameters have an insignificant effect on it. The results of the analysis of variance confirm this entry too.
Figure 8
The effect of flow rate variation with irreducible water saturation on the amount of recovered gas volume.
Figure 9
The effect of flow rate variation with permeability on the amount of recovered gas volume.
Figure 10
The effect of flow rate variation with pressure on the amount of recovered gas volume.
4.4. Interaction figure
Simultaneously studying the relationship between three parameters is a valuable way to analyze the effect of the parameters on a process. Figure 11 shows the interaction of two parameters, namely pressure and the flow rate of the injected gas, on the stored gas volume. It can be seen that by changing the flow rate of the injected gas, the distance between the two lines does not change much, and they remain parallel. This indicates that these two parameters have the minimum effect to each other on the stored gas volume[M.N.5] [R6] . It means that changing the flow rate of injected gas at different pressures has a same trend of effect on the stored gas volume. In addition, further investigation showed that there is not any significant interaction between permeability and irreducible water saturation.
Figure 11
Interaction of the flow rate of the injected gas and pressure on the stored gas volume.
4.5. Determination of optimum conditions and confirmatory test
The ability to optimize processes is one of the best advantages of the design of experiment methodologies. In this study, to be able to have the maximum storage volume efficiency in the natural gas storage process, all of the parameters affecting this process must be at an optimum level. Therefore, three different modes were considered in order to optimize the stored gas volume and the recovered gas volume. The first mode is maximizing the stored gas volume regardless of the recovered gas volume. The second mode is maximizing the recovered gas volume regardless of the stored gas volume; the third mode is simultaneously maximizing the stored gas volume and the recovered gas volume.
To this end, according to the selected ranges of each parameter, the best combination of factors is recommended as the optimum conditions by the model and is tabulated in Table 12.
To check the validity of the optimized conditions obtained by the model, the value of stored gas volume and the recovered gas volume obtained by the simulation were compared to the values predicted by the model of Design Expert software. A confirmation test was carried out by applying the optimum level of each factor. The predicted range of each response by considering a 95% confidence interval is [7.63111.9] for the stored gas volume in mode 1, [48.9664.46] for the recovered gas volume in mode 2, and [7.63111.95] and [48.9464.46] for the stored gas volume and the recovered gas volume in mode 3 respectively. The obtained results of running the optimal test are 8.6478 for the stored gas volume in mode 1, 52.75 for the recovered gas volume in mode 2, and 9.1068 for the stored gas volume and 54.99 for the recovered gas volume in mode 3. According to the results of the optimal test, it can be stated that the results of the verification test in each mode are close to its corresponding predicted value obtained by the statistical model, which confirms the validity of the obtained optimum test.
Table 12
Optimal conditions to obtain the highest response of the system in each mode.
Optimum value 
Parameters 
Mode 
0.85 
Flow rate (cc/min) 
1 
0.46 
S_{wirr} 

2044.12 
Pressure (psi)) 

1.75 
Permeability (md) 

0.82 
Flow rate (cc/min) 
2 
0.47 
S_{wirr} 

2135.61 
Pressure (psi)) 

1.09 
Permeability (md) 

0.9 
Flow rate (cc/min) 
3 
0.46 
S_{wirr} 

2254.89 
Pressure (psi)) 

1.54 
Permeability (md) 
4.6. Retention time
Gas was injected into the core and retained for four different times. This was conducted to investigate the effect of retention time on storage process. Gas was remained in the core for 0, 0.5, 2, and 10 hours, and it was then produced at a constant pressure. The results showed that as the retention time increases, the amount of the produced gas falls. The recovery of gas was between 46.31% and 55.76%., which shows the importance of retention time in the gas storage process. The reason for less gas production is gas dispersion inside the porous medium. Gas has high mobility and tends to move immediately after injection. In addition, water displacement and porous medium heterogeneity lead to gas trapping in some areas of the medium. Table 13 illustrates the results of the retention time tests.
Table 13
Effect of retention time on gas storage.
Recovery (%) 
Produced water (cc) 
Produced gas (cc) 
Time (hrs) 
55.76 
4.6354 
122.746 
0 
53.40 
4.6845 
116.920 
0.5 
46.76 
4.7092 
102.329 
2 
46.31 
4.7175 
101.334 
10 
5. Conclusions
In this study, the simulation of underground storage of natural gas in an aquifer and the investigation of parameters affecting this process are performed by Eclipse 300 software. The sensitivity analysis and grid independency test of numerical simulations were checked. The numerical results were compared with the experimental data, and good agreement was seen. Therefore, the uncertainty analysis of the numerical results was checked. At first, all of the parameters impacting on the process were determined; then, by using PlackettBurmanmethod, the most important factors were selected to plan a systematic study. Finally, Design Expert software with the help of RSM methodology was used to study the parameters influencing the process and to find the best combination of the factors which leads to the maximum stored and recovered gas volume. In overall, the following results were obtained:
Acknowledgments
The authors would like to thank Iran Natural Gas Storage Company (NGSC) for the financial support and StatEase, Minneapolis, USA for the provision of the Design Expert package.
Nomenclature
CV 
Coefficient of variance 
DOF 
Degree of freedom 
G 
Gravity 
K 
Permeability 
The standard average of the results at level j 

The total standard average 

Mw 
Molecular weight 
N 
The number of tests in PlackettBurman method 
N 
Number of levels 
P 
Pressure 
P_{c} 
Capillary pressure 
PV 
Pore volume of the injected fluid 
q 
Flow rate 
S 
Saturation 
Sum of squares 

S_{Wirr} 
Irreducible water saturation 
T 
Time 
Variance 

y+ 
Response of the system at the maximum level 
y 
Response of the system at the minimum level 
Z 
Height 
Greek letters 

Viscosity 

Density 

Porosity 

Subscripts 

w 
Water 
g 
Gas 