Optimizing Underground Natural Gas Storage Capacity through Numerical Modeling and Strategic Well Placement

Abstract

This study focuses on optimizing the storage capacity of an underground natural gas storage facility through numerical modeling and simulation techniques. The reservoir, characterized by an elongated dome structure, was discretized into approximately 16,000 cells. Simulations were conducted using key parameters such as permeability (10–70 mD) and porosity (12–26%) to assess the dynamics of gas injection and pressure distribution. The model incorporated core and petrophysical data to accurately represent the reservoir’s behavior. By integrating new wells in areas with storage deficits, the model demonstrated improvements in storage efficiency and pressure uniformity. The introduction of additional wells led to a significant increase in storage volume from 380 to 512 million Sm³ and optimized the injection process by reducing the storage period by 25%. The study concludes that reservoir performance can be enhanced with targeted well placement and customized flow rates, resulting in both increased storage capacity and economic benefits.

1. Introduction

The pressing global energy challenges and the need for sustainable human development demand the integrated development and utilization of energy resources, particularly oil and gas wells. Thus, conducting numerical studies on the factors influencing the stability of these wells is of critical importance, as they remain among the most essential energy sources [1].

Natural gas, with its abundance, versatility, and clean-burning properties, is playing an increasingly vital role in meeting the world’s energy demands [2,3]. According to forecasts from the International Energy Agency and the U.S. Energy Information Administration, the global demand for natural gas is expected to grow at an average annual rate of over 7.5% over the next 20 years, with industrialized nations seeing a growth rate exceeding 3% [4]. As time progresses, the demand for natural gas will continue to rise [5]. The steady energy supply from the pipeline system is effectively balanced with fluctuating market demands through the use of underground storage [6]. These storage reservoirs serve as gas supply warehouses, capable of meeting peak demand periods [7,8]. The underground natural gas storage market size was valued at USD 415.29 billion in 2023. The underground natural gas storage market industry is projected to grow from USD 430.20 billion in 2024 to USD 550.68 billion by 2032, exhibiting a compound annual growth rate (CAGR) of 3.13% during the forecast period (2024–2032) [9]. Natural gas storage is a complex process that must take into account the geology of the depleted deposit, the injection/extractive well equipment, and the actual storage process, which requires a large energy consumption. The main purpose of natural gas storage is to ensure the necessary consumption in the cold periods of the year, injecting the largest possible quantities of gas in the warm periods of the year, but let us not forget energy security, which has been debated more and more in recent years due to geopolitical events.

Natural gas is stored in underground geological formations, often repurposing former production fields into storage facilities (depleted gas reservoirs).

In some instances, storage sites are created in aquifer reservoirs that never previously held gas or within salt cavern reservoirs [10,11].

According to the International Gas Union, there are 692 underground natural gas storage (UGS) facilities operating globally. Of these, the majority (474) are situated in depleted natural gas fields. Additionally, 102 UGS facilities are located in salt caverns, while 41 are in depleted oil fields, and only 3 are in rock caverns [12].

The first gas storage facility at Dazhangtuo in the Dagang oil field, northern China, was established in 2004. Following this, five additional facilities (Ban-876, Banzhongnan, Banzhongbei, Ban-808, and Ban-828) were constructed by the end of 2006, with respective working gas volumes of 600, 145, 215, 612, and 250 million cubic meters (Vn). In Europe, several underground natural gas storage sites include Urziceni, Balaceanca, and Bilciuresti in Romania, with storage capacities of 350, 50, and 1310 million cubic meters (Vn), respectively. Other prominent storage sites include Bordolano in Italy (995 million cubic meters), Bernburg in Germany (971 million cubic meters), Bilche-Volytsko-Uherske in Ukraine (17,050 million cubic meters), Rough in the United Kingdom (17,050 million cubic meters), and Kasimovskoe in Russia (10,250 million cubic meters). In North America, notable projects include Fink-Kennedy-Lost Creek, with a capacity of 2535 million cubic meters (Vn) [13].

The efficiency of a gas storage facility, in terms of withdrawal and injection capacity, depends on key technical factors such as reservoir pressure, the physical characteristics of the geological formation, and the number of wells. Depleted gas deposits are complex structures located on large areas that are exploited by means of numerous injection/extraction wells and which are essential in optimizing the storage process. The efficiency of this process depends on the permanent knowledge of the dynamics of the processes and the distribution of the gas cushion. The maximum gas injection capacity is typically maintained throughout the entire storage cycle, whereas the maximum withdrawal capacity is only sustainable up to a certain reservoir pressure, which gradually decreases as gas is depleted [14]. A study conducted by Moradi used the commercial simulator Eclipse 300 to model the impact of natural gas injection and production on gas relative permeability and condensate saturation in both the near-wellbore region and areas farther from the wellbore. The study revealed that injecting a larger volume of gas in the initial injection phase could effectively prevent the formation of condensate near the wellbore area [15]. Shin and Lee [16] investigated the effectiveness of an operational scenario for an offshore gas condensate storage reservoir composed of five layers using the Eclipse 300 compositional simulator. This model enabled a detailed analysis of pressure variations within the reservoir, confirming a significant pressure difference influenced by phase behavior during operations. By closely examining the compositional changes in the gas, the simulation provided valuable technical data that supported the selection of optimal operation capacity and methods.

The design of gas storage operations hinges on the reservoir’s location and performance [17]. Decisions aimed at maximizing gas volume and deliverability, such as determining where to drill new wells and identifying which wells may need to be abandoned, must be grounded in a thorough analysis of the reservoir’s conditions, supported by accurate simulations of its behavior.

Failure to capitalize on the potential of gas deposits can have dramatic consequences for the increase in prices and the decrease in Romanian competitiveness.

The location of large consumers at the ends of the systems, far from the main gas sources, represents a significant disadvantage because the reaction time of the extraction–transportation system is very high compared to the consumption needs [18]. During the cold season, these capacities become insufficient, especially at peak consumption times.

This approach ensures that operations and investments are executed effectively [19]. Neglecting these critical parameters can lead to significant risks, including gas leakage and operational failure. Such leakage may occur due to factors like gas migration, entrapment in inaccessible reservoir regions, or the fingering phenomenon in the displacing gas phase [20].

Mazarei et al. [21] probed and optimized the process of gas storage operations in a depleted fractured gas condensate reservoir in Iran. Their findings concluded that determining the optimal volume of working gas not only maximizes gas and condensate production but also ensures that the planned production rate for the recovery season is achievable. Additionally, they found that the closer the composition of the injected gas is to that of the pipeline gas, the better the mixing with the reservoir gas, leading to higher gas condensate recovery. Their study also highlighted that adding more wells significantly increased cumulative gas and condensate production rates while reducing cumulative water production, thereby enhancing overall recovery efficiency.

Most of the natural gas deposits located on the territory of Romania are depleted, which means that the favorable conditions have been met to develop underground gas deposits that can ensure the necessary internal consumption but which can also lead to the development of favorable external markets from an economic point of view, Romania potentially becoming a strategic exporter of natural gas, being favorable in the current geopolitical context.

When storing, one of the most important technical and economic factors of deposits is represented by the ratio between gas cushion and useful capacity [22]. Through numerical modeling, we created the computer model that reproduces the geometry of the deposit according to the properties of the rocks, and thus we were able to define the movement of fluids in the deposit.

Numerical simulations are carried out on the computer model of the deposit with the aim of predicting the exploitation of the deposit over different time intervals. Through these procedures, the phases of the gas storage process in depleted deposits can be simulated.

In the current paper, the authors describe a numerical nonstationary two-dimensional model that would allow the simulation of the storage processes by wells in time and space, thus obtaining a dynamic spatial image of the gas distribution in storage. Gas2Storage v1 software can be used to check up the nominations depending on the possibilities of the storage. The optimization of the storage process along with the redistribution of gas in the reservoir based on the possibility of using new wells will ensure the best decisions for management.

2. Methodology

In order to carry out this study, data were needed regarding the average parameters of the reservoir, the initial pressure, and the volume of the existing gas cushion in the reservoir, as well as the average seasonal injection rate and the average storage and equilibrium period of the reservoir. Leveraging this dataset received from the storage company, simulations were conducted and analyzed to assess the introduction of natural gas in the reservoir in optimal conditions and with maximum efficiency.

The working methodology engaged to write this article is visually depicted in Figure 1. The steps taken are:

• Choosing the optimization of the storage process to be enhanced based on the new findings in the literature and the development of the simulation software.
• Modeling of the reservoir topology and injection wells with variable quantities of natural gas.
• Constructing working scenarios for new well integration involving the considering various factors, including the average seasonal volumes and reservoir distribution.
• Performing a results analysis of the considered scenarios emphasizing the dynamics of the stored gas and distribution in the reservoir in the scope of diminishing the energy used, as well as an increase in volumes.
• Holding a comprehensive discussion of the results and presenting their implications and significance for reservoir capacity enhancement.
• Outlining the potential for future research in the conclusions of the study based on its findings and insights.

The underground natural gas storage facility has a capacity of 300 million Sm³ with an initial and constant reservoir pressure of 30 bar. The reservoir model, which is based on the geological structure’s engineering, is discretized into approximately 16,000 cells. Structurally, the reservoir takes the form of an elongated dome oriented in the NW-SE direction, with consistent depth and relatively low inclinations. It is treated as a single hydrodynamic unit, and due to its unique lithological properties, a vertical saturation limit has been established along the −900 m isobath.

Core analysis shows porosity values between 20% and 26%, aligning with estimates from regional porosity maps. However, petrophysical analysis indicated lower porosity values, ranging from 12% to 18%, which are significantly below the core data results. Permeability averages between 10 and 70 mD, and for simulation purposes, an average permeability of 10 mD and a porosity of 20% were assumed. As a result, pressure increases are more pronounced near the wells and gradually dissipate towards the reservoir’s boundaries.

3. Research on Increasing Gas Reservoir Storage

For the precise simulation of the dynamic parameters and the distribution in the reservoir, a numerical model was developed, being of maximum interest for analyzing the pressure matrix and gas volume distribution.

3.1. Presentation of the Simulation Model

The two-dimensional movement of gasses in a homogeneous porous medium relative to pressure is defined by Equation (1).

$${\displaystyle \frac{\partial }{\partial x}}\left({\displaystyle \frac{k·h}{\mu ·z}}·{\displaystyle \frac{\partial p^2}{\partial x}}\right)+{\displaystyle \frac{\partial }{\partial y}}\left({\displaystyle \frac{k·h}{\mu ·z}}·{\displaystyle \frac{\partial p^2}{\partial y}}\right)=2·m·h·{\displaystyle \frac{\partial }{\partial \tau }}\left({\displaystyle \frac{p}{z}}\right)$$

(1)

The field of integration is the deposit (Figure 2).

In order for Equation (1) to be solved numerically, it is necessary to technically define the domain of integration. Thus, the deposit domain is divided into blocks of size $dx\times dy\times h$.

A new variable is introduced, called pseudo pressure:

$$u={\displaystyle \frac{p^2}{\mu ·z}}$$

(2)

with the partial derivatives obtained:

$$\begin{gathered} {\displaystyle \frac{\partial u}{\partial x}}={\displaystyle \frac{1}{\mu ·z}}·{\displaystyle \frac{\partial p^2}{\partial x}} \\ {\displaystyle \frac{\partial u}{\partial y}}={\displaystyle \frac{1}{\mu ·z}}·{\displaystyle \frac{\partial p^2}{\partial y}} \\ {\displaystyle \frac{\partial u}{\partial \tau }}={\displaystyle \frac{1}{\mu ·z}}·{\displaystyle \frac{\partial p^2}{\partial \tau }}={\displaystyle \frac{2·p}{\mu ·z}}·{\displaystyle \frac{\partial p}{\partial \tau }} \end{gathered}$$

(3)

For left members of Equation (1), using (3), obtained:

$${\displaystyle \frac{\partial }{\partial x}}\left({\displaystyle \frac{k·h}{\mu ·z}}·{\displaystyle \frac{\partial p^2}{\partial x}}\right)+{\displaystyle \frac{\partial }{\partial y}}\left({\displaystyle \frac{k·h}{\mu ·z}}·{\displaystyle \frac{\partial p^2}{\partial y}}\right)={\displaystyle \frac{\partial }{\partial x}}\left(k·h·{\displaystyle \frac{1}{\mu z}}·{\displaystyle \frac{\partial p^2}{\partial x}}\right)+{\displaystyle \frac{\partial }{\partial y}}\left(k·h·{\displaystyle \frac{1}{\mu z}}·{\displaystyle \frac{\partial p^2}{\partial y}}\right)$$

(4)

For $k$ and $h$ with median values, Equation (4) becomes:

$$k·h·{\displaystyle \frac{\partial }{\partial x}}\left({\displaystyle \frac{\partial u}{\partial x}}\right)+k·h·{\displaystyle \frac{\partial }{\partial y}}\left({\displaystyle \frac{\partial u}{\partial y}}\right)=k·h·\left({\displaystyle \frac{{\partial }^{2}u}{\partial x^2}}+{\displaystyle \frac{{\partial }^{2}u}{\partial y^2}}\right)$$

(5)

For the right members of Equation (1), using (3), obtained:

$${\displaystyle \frac{\partial }{\partial \tau }}\left({\displaystyle \frac{p}{z}}\right)={\displaystyle \frac{1}{z}}·{\displaystyle \frac{\partial p}{\partial \tau }}+p·{\displaystyle \frac{\partial z}{\partial \tau }}={\displaystyle \frac{1}{z}}·{\displaystyle \frac{\partial p}{\partial \tau }}-p·{\displaystyle \frac{\partial z}{\partial p}}·{\displaystyle \frac{\partial p}{\partial \tau }}=\left({\displaystyle \frac{1}{z}}-p·{\displaystyle \frac{\partial z}{\partial p}}\right)·{\displaystyle \frac{\partial p}{\partial \tau }}={\displaystyle \frac{p}{z}}·\left({\displaystyle \frac{1}{p}}-{\displaystyle \frac{1}{z}}·{\displaystyle \frac{\partial z}{\partial p}}\right)·{\displaystyle \frac{\partial p}{\partial \tau }}$$

(6)

Using notation for compressibility:

$$\beta ={\displaystyle \frac{1}{\rho }}{\displaystyle \frac{\partial \rho }{\partial p}}={\displaystyle \frac{z}{p}}{\displaystyle \frac{\partial }{\partial p}}\left({\displaystyle \frac{\rho }{z}}\right)={\displaystyle \frac{1}{p}}-{\displaystyle \frac{1}{z}}{\displaystyle \frac{\partial z}{\partial p}}$$

(7)

Equation (6) becomes:

$$2·m·h·{\displaystyle \frac{\partial }{\partial \tau }}\left({\displaystyle \frac{p}{z}}\right)=m·h·\beta ·\mu ·\left({\displaystyle \frac{2·p}{\mu ·z}}·{\displaystyle \frac{\partial p}{\partial \tau }}\right)=m·h·\beta ·\mu ·{\displaystyle \frac{\partial u}{\partial \tau }}$$

(8)

With this notation, considering the division of the domain, Equation (1) written in finite differences becomes:

$$k·\left({\displaystyle \frac{\partial u}{\partial x^2}}+{\displaystyle \frac{{\partial }^{2}u}{\partial y^2}}\right)=m·\mu ·\beta ·{\displaystyle \frac{\partial u}{\partial \tau }}$$

(9)

Using explicit methods for differential equations:

$$\begin{gathered} {\displaystyle \frac{{\partial }^{2}u}{\partial x^2}}\cong {\displaystyle \frac{u_{i+1,j}^n-2u_{i,j}^n+u_{i-1,j}^n}{\Delta x^2}} \\ {\displaystyle \frac{{\partial }^{2}u}{\partial y^2}}\cong {\displaystyle \frac{u_{i,j+1}^n-2u_{i,j}^n+u_{i,j-1}^n}{\Delta y^2}} \\ {\displaystyle \frac{\partial u}{\partial \tau }}\cong {\displaystyle \frac{u_{i,j}^{n+1}-u_{i,j}^n}{\Delta \tau }} \end{gathered}$$

(10)

obtained:

$$k·\left({\displaystyle \frac{u_{i+1,j}^n-2u_{i,j}^n+u_{i-1,j}^n}{\Delta x^2}}+{\displaystyle \frac{u_{i,j+1}^n-2u_{i,j}^n+u_{i,j-1}^n}{\Delta y^2}}\right)=m·\beta ·\mu ·\left({\displaystyle \frac{u_{i,j}^{n+1}-u_{i,j}^n}{\Delta \tau }}\right)$$

(11)

Using notations:

$$\begin{gathered} {\alpha }_x={\displaystyle \frac{k·\Delta \tau }{m·\beta ·\mu ·\Delta x^2}} \\ {\alpha }_y={\displaystyle \frac{k·\Delta \tau }{m·\beta ·\mu ·\Delta y^2}} \end{gathered}$$

(12)

Equation (11) with the notes above becomes:

$$k·\left({\displaystyle \frac{u_{i+1,j}^n-2u_{i,j}^n+u_{i-1,j}^n}{\Delta x^2}}+{\displaystyle \frac{u_{i,j+1}^n-2u_{i,j}^n+u_{i,j-1}^n}{\Delta y^2}}\right)=m·\beta ·\mu ·\left({\displaystyle \frac{u_{i,j}^{n+1}-u_{i,j}^n}{\Delta \tau }}\right)$$

(13)

After the processing of Equation (13):

$$u_{i,j}^{n+1}={\alpha }_x·\left(u_{i+1,j}^n+u_{i-1,j}^n\right)+{\alpha }_y·\left(u_{i,j+1}^n+u_{i,j-1}^n\right)+u_{i,j}^n·\left(1-2·{\alpha }_x-2·{\alpha }_y\right)$$

(14)

Since the variables at time $“n+1”$ are not known, Equation (14) generates a system of equations from which the variable $u_{i,j}$ is determined in all nodes of the network at any time.

3.1.1. Initial Conditions

These conditions represent the pressure values in the network nodes at the beginning of the integration, basically the initial pressure in the reservoir.

$$u_{i,j}={\displaystyle \frac{p_{i,j}^2}{\mu ·z}}={\displaystyle \frac{p_{init}^2}{\mu ·z}}$$

(15)

3.1.2. Boundary Conditions

In the deposit, the integration blocks can be of several kinds (Figure 3).

The pressure required for the gas injection is superimposed on the existing pressure in the deposit, which we call the pressure at the base of the well and denote it by $p_{bw}$.

The pressure generated by the massive flow of injected gas, accumulated with the pressure at the base of the well ($p_{bw}$), is:

$$p_{ij}=p_{bw}+\mathsf{\varphi }·{\displaystyle \frac{\mu }{2·\pi ·k·h·\rho }}$$

(16)

where:

$$p_{bw}={\displaystyle \frac{\left(p_{i+1,j}+p_{i-1,j}+p_{i,j+1}+p_{i,j-1}\right)}{4}}$$

(17)

In each block where a well is placed, the pressure condition resulting from the injected flow rate is applied (Figure 3).

3.1.3. Domain Boundary Conditions

Given that we are dealing with a gas field, its boundary is considered impermeable if block $[i,j]$ is a boundary block (Figure 4). The gas velocity between it and the adjacent block in the reservoir must be zero (18).

$$\begin{gathered} w_x=-{\displaystyle \frac{k}{\mu }}·{\displaystyle \frac{\partial p}{\partial x}}\cong -{\displaystyle \frac{k}{\mu }}·{\displaystyle \frac{{\left(p_{i,j}^{n+1}\right)}_f-p_{i+1,j}^{n+1}}{\Delta x}}=0 \\ {w}_y=-{\displaystyle \frac{k}{\mu }}·{\displaystyle \frac{\partial p}{\partial y}}\cong -{\displaystyle \frac{k}{\mu }}·{\displaystyle \frac{{\left(p_{i,j}^{n+1}\right)}_f-p_{i,j+1}^{n+1}}{\Delta y}}=0 \end{gathered}$$

(18)

Using (18) obtained:

$$\begin{gathered} {\left(p_{i,j}^{n+1}\right)}_f=p_{i\pm 1,j}^{n+1} \mathrm{on} \mathrm{X} \mathrm{direction} \\ {\left(p_{i,j}^{n+1}\right)}_f=p_{i,j\pm 1}^{n+1} \mathrm{on} \mathrm{Y} \mathrm{direction} \end{gathered}$$

(19)

where the index $f$ refers to the boundary.

The integration process is carried out for each time step as follows:

• The initial conditions are placed on all reservoir blocks.
• The boundary conditions are set.
• The system of equations is generated according to (14).
• The system of equations is solved, resulting in the gas pressure in each block.
• Processes 2, 3, and 4 are repeated until the sting period is over.

3.1.4. Calculation of the Gas Mixture Process on the Deposit

If the injected gasses have a different composition, then the model must take this into account. For this purpose, we define in each block the gas composition (${CMP}_{i,j}$), which can be molar or mass. Mass composition is useful because mass balance is used in the calculation. The composition calculation is performed simultaneously with the reservoir pressure calculation, usually for the same time step. The composition of block $[i,j]$, $M_{i,j}^{n+1}\cdot g\left(k\right)$ is considered because it does not change much for a calculation step $∆\tau$.

After the reservoir pressure calculation, the mass balance is made for all the blocks in order to determine the new compositions of the block.

Knowing that the mass flow rate has the expression:

$$\dot{m}=\mathsf{\rho }·\dot{V}=\mathsf{\rho }·w·A$$

(20)

using the notations in Figure 5, we can define the mass fluxes as:

$$\begin{gathered} \dot{m_1}=-{\displaystyle \frac{k}{\mu }}·\left({\displaystyle \frac{p_{i,j}-p_{i,j+1}}{{∆}_y}}\right)·∆x·h·{\rho }_{i,j+1} \\ \dot{m_2}=-{\displaystyle \frac{k}{\mu }}·\left({\displaystyle \frac{p_{i,j}-p_{i+1,j}}{{∆}_x}}\right)·∆y·h·{\rho }_{i+1,j} \\ \dot{m_3}=-{\displaystyle \frac{k}{\mu }}·\left({\displaystyle \frac{p_{i,j}-p_{i,j-1}}{{∆}_Y}}\right)·∆x·h·{\rho }_{i,j-1} \\ \dot{m_4}=-{\displaystyle \frac{k}{\mu }}·\left({\displaystyle \frac{p_{i,j}-p_{i-1,j}}{{∆}_x}}\right)·∆y·h·{\rho }_{i-1,j} \end{gathered}$$

(21)

The mass balance on the well is expressed as follows:

$$\dot{m}=\dot{m_1}+\dot{m_2}+\dot{m_3}+\dot{m_4}$$

(22)

The gas mass of block $[i,j]$ is determined with the value of the density at the previous time step:

$$m_b=∆x·∆y·h·{\rho }_{i,j}$$

(23)

Using the mass balance of the block, its new composition can be calculated for every $k$ component as follows:

$$M_{i,j}^{n+1}\cdot g\left(k\right)={\displaystyle \frac{\sum _{\begin{array}{c}s=1\\ m_s>0\end{array}}^4{m}_s\cdot M_{i,j}^n\cdot g\left(k\right)+m_b \cdot M_{i,j}^n\cdot g\left(k\right)}{m_b+\sum _{\begin{array}{c}s=1\\ m_s>0\end{array}}^4{m}_s}}$$

(24)

Only positive mass fluxes (entering the block) are considered because they bring the new composition. Negative gas flows leave the block with its composition and will be taken into account in adjacent blocks.

3.2. Presentation of the Gas2Storage Simulation Software

Building the Gas2Storage simulation software started from the real geometry of a deposit with a flat shape, and because it has a relatively constant thickness, a two-dimensional model was created. Figure 6 shows an image of the deposit contour loading and the location of available and future wells. These probes can be used in whole or in part for the gas storage phase.

The definition of the computer map of the domain is performed by assigning different attributes for the points inside the deposit, the border, the position of the probes, and the points outside the deposit. They help to carry out the integration process and introduce the conditions on the contour and in the well, allowing to organize the calculations for integration domains with different geometries.

The gas composition used is from actual chromatographic bulletins given by the operator associated with the entry point from the transmission system. The gas used in the simulation has 92.4 Methane, 2.1 Ethane, 1.5 Propane, 0.9 N—Butane, 0.5 ISO—Butane, 0.3 Pentane, 0.5 Hexane, 0.6 Nitrogen, and 1.2 Carbon dioxide.

To define the various scenarios applied to the simulation process, we need to define initial conditions and boundary conditions on the boundary and in the wells. The initial conditions reflect the state from which the simulation starts. The amount of gas initially existing in the reservoir has a pressure that we can call the initial pressure. The boundaries are considered impermeable in the case of gas deposits, with the gas velocity at the boundary points being zero. The boundary condition in the well is represented by the injected or extracted flow rate. This depends on the dimensions of the well, the reservoir pressure, and the permeability of the formation. In the scenario, flow through the well is imposed, and the model will result in the injection/extraction pressure that ensures the circulation of that flow through the well. The coincidence of the loading curves of the deposit and the model shows that the numerical simulator is calibrated. The model is able to present a spatial image of gas pressure distribution in the reservoir. The peaks on the surface represent the pressures in the wells (Figure 7).

During the quiescence phase, the gas migrates towards the boundary of the field, but due to the location of the wells, a significant portion of the field remains unused. After this phase, it is possible to observe the uniformity of the pressures in the area of the injection wells, and new wells can be optimally placed and distributed where there are deficient areas in order to increase the storage capacity of the deposit.

In order to increase the benefits obtained within the gas storage process, the amount of gas that can be stored in the respective deposit must be maximized by digging new wells and by determining a value of the injected flow rate.

The analysis of the spatial distribution of the reservoir gasses at the end of the injection period highlights an important non-uniformity. In order for the filling of the deposit to be uniform, due to the low filtration speed, a quiet period in the order of months would be needed, much too long for the storage process.

As the gas storage takes place over a limited period of time, the gasses move in the reservoir for small distances around the wells, leaving portions of the reservoir unfilled. To correct this problem in the simulation, new probes were introduced located in the area with a filling deficit, the green area in Figure 8.

3.3. Outcomes of the Simulation

The results of the numerical modeling highlighted the areas that have a storage deficit, thus making it possible to determine the location of the new wells. The gas flows injected through all wells in the deposit had equal average values (4166.67 Sm³/h) in order to highlight the results obtained after the process (Figure 9).

The presentation at various points in time gives us clear indications of the increase in the amount of stored gas. In the following, the comparative results obtained from every 30 days of the entire storage process are presented (Figure 10, Figure 11, Figure 12 and Figure 13). The every 30 days means that we have performed simulations and present the results every 30 days to see the differences in time. The simulation results for the actual situation show a parabolic increase in average pressure in the wells and also on the average pressure in the reservoir during time starting from 47 bar and 33 bar as is shown in Figure 9-left for 30 days of injection, 52 bar and 34 bar after 60 days (Figure 10-left), 57 bar and 39 bar after 90 days (Figure 11-left), 42 bar and 20.2 bar after 90 days (Figure 12-left), and finally 61 bar and 42 bar after 120 days (Figure 13). As it was expected, the flow through the pore is limited because of the reservoir parameters. If we analyze the right part of Figure 10, Figure 11, Figure 12 and Figure 13 where we present the values with the new wells compared with the right part of the figures, we can observe that for each 30 days the average pressure in the wells decreases with 1 bar, and also the average pressure in the reservoir increases with 1 bar, ultimately obtaining a difference of 6 bar between them.

The pressure increases faster in the central area due to the high density of the probes and the low gas filtration rate. The contribution of the probes located in the newly introduced deficit regions in the modeling is observed.

The gas distribution in the reservoir at the end of the injection shows (Figure 14a) a more uniform loading of the reservoir compared to the original situation (Figure 14b).

If we analyze this distribution compared to the case where we did not have new wells and after the end of the quiet period (Figure 15), the same clear difference in distribution is obtained, which leads to an increase in the average pressure on the field and, of course, to the storage of a larger amount of gas.

The analysis of the cumulative data of the storage processes, Figure 16, shows that the results of the introduction of new wells are also felt in the values obtained; thus, the injected flow rate is approximately 512 million Sm³ compared to 380 million Sm³ in the initial version, and the average pressure on the deposit has increased although the average pressure of the wells decreases.

4. Discussions

In general, the most economical way to increase the delivery capacity of the deposits is to increase the volume of the gas cushion, which automatically implies an increase in the pressure level, generating an increased energy support on the basis of which larger quantities can be obtained over the same period of time extraction.

The numerical simulator was calibrated with SCADA data for both the injection and extraction phases. It should be noted that this process must be carried out periodically.

For the case where the reservoir is in the quiet phase, i.e., no gas is injected, no gas is extracted, the pressure in the reservoir tends to equalize.

To calibrate the Gas2Storage software, production data for 3 months of injection as well as the initial reservoir pressure were used. A comparison was made between the cumulative injected pressure in the model and the one realized in the real deposit.

Since the gas consumption process is not uniform and depends a lot on temperatures and industry, the possibilities of requesting the deposit by temporarily increasing the injected flows in a certain period must be known.

The research is carried forward by making some curves that represent the connection between the injected flow rates, the amount of gas, and the average pressure in the reservoir. This allows the quick determination of the average daily flow per well that must be injected in order to obtain a certain pressure in the deposit at the end of the cycle.

For the most optimal loading of the gas deposit, the spatial distribution of the gasses in the deposit must be known—something that can be visualized and calculated by analyzing the pressure field. The more uniform filling we obtain, the more we increase the storage performance. These things were also observed by Gligor et al. [23]. The main cause of the uneven distribution of gas in the deposit is the location of the wells, dug according to the needs of the extraction process, and the introduction of new wells in the deficient areas can lead to significant increases in the stored gasses. As the wells were drilled a long time ago, the optimization process was not so evolved as presented in more recent articles, one of which was our work Production Forecasting at Natural Gas Wells [24].

As can be seen from the presented figures, the introduction of new wells in the underprivileged areas of the deposit has led to an increase in both the amount of stored gas and the average pressure on the deposit due to the much more efficient distribution of gasses on the surface of the deposit. This was achieved even with lower injection pressures, which reduces compression costs and therefore increases the efficiency of the process.

Another advantage of the redistribution of the injected flows would be the reduction in the storage period by 25%, which would help in the context of discussions related to nominations and gas price variability.

Future analyses will target a much more accurate distribution of injected volumes by using customized flow rates on each well, depending on the position in the reservoir and the distribution and migration of gasses within it. In addition to verifying the increase in gas storage capacity based on the number of new wells, this can also be provide an economic calculation.

5. Conclusions

The operating mode of the wells to fulfill the imposed requirements was defined based on numerical simulations that respected the reservoir model and the nominations.

An original two-dimensional, non-stationary, variable-geometry numerical model for a gas reservoir was created. It allows the daily determination of the matrix of pressures, gas saturation, and gas volumes at any point in the reservoir and its representation in 3D graphic form.

The developed Gas2Storage software can be used for verifying and accepting nominations, finding solutions for crisis situations, or verifying depleted deposits proposed to be transformed into underground storage deposits.

The analysis of the obtained results leads to the veracity of the storage optimization process because the average pressure on the deposit has increased although the average pressure of the wells is decreasing, starting at 47 bar in the wells and 33 bar in the reservoir after 30 days of injection, rising to 52 bar and 34 bar after 60 days, and finally reaching 61 bar in the wells and 42 bar in the reservoir after 120 days, which means a much lower energy consumption per well in conjunction with the increase in the amount of gas stored in the deposit.

The increase in the number of wells, although it requires large investments, if it is performed correctly can lead to additional income due to the much larger volumes of natural gas stored, but also to the provision of storage/extraction possibilities from various areas of the deposit depending on their distribution, in favorable periods financially, depending on gas prices and needs.

As the next step, a technical and economic study should be undertaken encompassing both the establishment of storage capacity and the amortization period associated with this technical objective.