A Novel Performance Evaluation Method for Gas Reservoir-Type Underground Natural Gas Storage

Abstract

The regulation of the seasonal energy supply for natural gas and the storage of fossil energy are important to society. To achieve it, storing a large amount of natural gas in porous underground media is one of the government’s choices. Due to the successful lesson learned from the oil and gas industry, natural gas storage in underground porous media has been regarded as the most potential long−term energy storage method. In this paper, we developed a new workflow to evaluate the performance of gas reservoir−type underground natural gas storage (UGS). The theoretical background of this workflow includes the correction of the average formation pressure (AFP) and gas deviation factor by error theory and the analytical mathematical model of UGS wells. The Laplace transform, line source function, and Stehfest numerical inversion methods were used to obtain pressure solutions for typical vertical and horizontal wells in UGS. The pressure superposition principle and weighting method of the gas injection−withdrawal rate were used to obtain the AFP. Through the correction of the AFP and gas deviation factor in the material balance equation, the parameters for inventory, effective inventory (the movable gas volume at standard condition), working gas volume (the movable gas volume is operated from the upper limit pressure to the lower limit pressure), and effective gas storage volume (the available gas storage volume at reservoir condition) were determined. Numerical data from the numerical simulator was used to verify the proposed model pressure solution. Actual data from China’s largest Hutubi UGS was used to illustrate the reliability of the proposed workflow in UGS performance evaluation. The results show that large−scale gas injection and withdrawal rates lead to composite heterogeneity in gas storage wells. The nine injection and production cycles’ pressure and effective inventory changes from Hutubi UGS can be divided into a period of rapid pressure rise and a period of slow pressure increase. The final AFP is 32.8 MPa. The final inventory of the Hutubi UGS is 100.1 × 10⁸ m³, with a capacity filling rate (the ratio of effective inventory to designed gas storage capacity) of 93.6%. The effective inventory is 95.3 × 10⁸ m³, and the inventory utilization ratio (the ratio of effective inventory to inventory) is 95.2%. The working gas volume is 40.3 × 10⁸ m³. This study provides a new method for inventory evaluation of the gas reservoir−type UGS.

1. Introduction

The Paris Climate Agreement calls for the world to reduce carbon emissions to nearly 95% of their 1990 level by 2050 [1,2]. The Chinese government has also proposed the goal of carbon neutrality, to achieve “zero emissions” of carbon dioxide by 2060 [3,4,5,6]. Therefore, the transition from traditional fossil energy represented by coal and oil to new energy has become imminent, and this transition has put forward higher requirements for new, large−scale energy storage technologies [7,8]. Starting from the storage medium, current energy storage technologies include mechanical energy storage, electrical energy storage, electrochemical energy storage, thermal energy storage, and chemical energy storage [9,10]. Compared with other energy storage technologies, mechanical energy storage technologies represented by pumped storage, compressed air energy storage, and underground gas storage have been favored by more and more countries and regions due to their excellent safety, economy, and long−term performance [11,12,13]. Additionally, in order to better balance the seasonality and regionality of energy supply and demand throughout the year, underground natural gas storage (UGS) has become the government planning task for various countries [14,15].

Target UGS regional structures include depleted oil and gas reservoirs, salt caverns, aquifers, and abandoned underground mines [16]. Due to the characteristics of geological structure and the traps sealing property, depleted oil and gas reservoirs are currently the mainstream way of using UGS [17,18]. As shown in Figure 1, 67% of the 689 UGS projects in operation around the world in 2018 came from depleted oil and gas reservoirs [19,20]. These 689 operating gas storage projects do not include EOR−related components. The Chinese government has also been committed to the construction of UGS in the past few decades, and the vast majority of them belong to this type of depleted oil and gas reservoirs. In 1970, Daqing Oilfield from China National Petroleum Corporation started the construction of China’s first UGS project [21,22]. The Dazhangtuo UGS was completed and put into operation in 2000. It was the first commercially operated UGS in China. [23]. Jintan UGS, as China’s first salt cavern−type UGS, began production in 2018 [24]. In order to meet the purpose of west−east natural gas transmission, China National Petroleum Corporation started the Hutubi UGS construction in 2013 and it was designed to be the largest UGS in China [25].

Similar to the SEC reserves of oil and gas reservoirs, the effective inventory and working gas volume of UGS are important parameters for evaluating the company’s performance and directly determine the capital value of the gas storage company [26]. The successful lesson from gas reservoir development shows that the material balance method, the numerical simulation method, and the volumetric method are common methods used to calculate the geological reserves of gas reservoirs [27,28,29,30]. This also provides a reference for the effective inventory and working gas volume calculation in the field of gas storage. Many scholars have used commercial numerical simulators to evaluate the gas storage capacity during UGS [31,32,33,34,35,36,37]. Although the numerical simulator can simulate the entire gas injection−withdrawal cycle of the gas storage, the complexity of the composition model and geological model brought by gas injection further increases the scene limitations of model application. In particular, the alternation of additional injection−withdrawal cycles in the gas storage brings more computational cost to history matching in the numerical simulation [38,39,40]. The volumetric method can directly calculate the theoretical gas storage volume and capacity from the geological perspective. However, this value is static, and it ignores the dynamic behaviors of gas injection and withdrawal cycles. The material balance method is a common method in gas reservoir engineering to calculate the oil and gas reservoirs geological reserves, and it has been used in the field of gas storage [41,42,43]. Through the pressure and deviation factor ratios in the original state and the average state, the oil and gas reservoir reserves can be obtained by the Y−intercept. Therefore, the prerequisite for the application of the material balance method is that the balanced average formation pressure (AFP) at the end of the injection and withdrawal stages is known. However, there is still a lack of methods to determine the balanced AFP in gas storage. The current material balance methods in UGS usually only use static pressure data at a certain time point.

To narrow this gap, this paper proposes a performance evaluation method for gas reservoir−type UGS. The mathematical theoretical basis includes the correction of the AFP and deviation factor by error theory and the analytical mathematical model of gas storage wells. The gas diffusion equation in UGS was solved by the Laplace transform method and the Stehfest numerical inversion method. The point source and line source function methods were used to obtain bottom−hole pressure solutions for vertical and horizontal wells in UGS. The reliability of the pressure solution was verified by data from commercial numerical simulation software. Furthermore, the proposed method was used to evaluate the inventory of Hutuobi UGS.

2. Methodology

In this section, the mathematical models for vertical and horizontal gas storage wells were respectively established to describe the flow behaviors in UGS. The reservoir domain consists of two regions with different flow capacities. The point source, line source function, and Laplace transform methods were used to obtain the pressure solutions of vertical and horizontal well models in Laplace space. The wellbore storage and the skin effects were considered in the method proposed by Mukherjee and Economides [44]. The Stehfest inversion method was used to obtain the bottomhole pressure solution in the real domain [45].

2.1. Physical Model

The physical model of vertical and horizontal wells in UGS is shown in Figure 2. The near−wellbore zone with high flow capacity was caused by high−speed gas injection and withdrawal rates, and the outer region refers to the initial formation. The UGS was considered to have a consistent formation thickness and a uniform initial pressure and temperature distribution. Horizontally, the physical model of UGS was infinite. The inner and outer regions were assumed to be continuous media. The vertical and horizontal wells were located in the center of the study area. The duration of the test well process is usually short compared to the long duration of the gas injection and withdrawal processes. Therefore, the inner region radius was assumed to be a constant value $r_f$. The interface between inner and outer boundaries satisfied the conditions for equal pressure and unequal flow rates. The remaining assumptions are as follows:

(1)

Gas in the reservoir was the single−phase, compressible fluid, and it obeyed Darcy’s law.

(2)

The vertical and horizontal wells in UGS were involved in the gas injection and withdrawal processes at a constant rate.

(3)

The effects of gravity and temperature on gas flow were ignored.

(4)

The UGS formation had a closed top−to−bottom interface in the vertical direction.

(5)

The effects of wellbore storage and skin−on−gas flow during gas injection and withdrawal were considered.

(6)

Due to the gas’s compressibility, the pseudo−pressure method proposed by Al−Hussainy was used to describe the natural gas flow [46].

$$m=2{\displaystyle {\int }_{p_0}^p\frac{p}{\mu Z}}\mathrm{d}p$$

(1)

where $m$ is the pseudo−pressure, $p$ is the pressure, $\mu$ is the gas viscosity, $Z$ is the gas deviation factor, and $p_0$ is the reference pressure. According to Boyle’s law, the deviation factor can be expressed as Equation (2) [47].

$$Z\left(p,T\right)=\frac{p}{{\rho }_{m}RT}$$

(2)

where ${\rho }_m$ is the molar density, $T$ is the reservoir temperature, and $R$ refers to the gas constant.

Figure 2. Physical model of UGS: (a) schematic diagram of gas reservoir−type underground natural gas storage (modified from Geostock [48]); (b) the horizontal and vertical wells with heterogeneous features in underground natural gas storage.

Figure 2. Physical model of UGS: (a) schematic diagram of gas reservoir−type underground natural gas storage (modified from Geostock [48]); (b) the horizontal and vertical wells with heterogeneous features in underground natural gas storage.

2.2. Mathematical Model

With the physical model and its assumptions, the governing equations for the vertical and horizontal wells in UGS can be given as:

$$\frac{{\partial }^2{m}_1}{\partial r^2}+\frac{1}{r}\frac{\partial m_1}{\partial r}=\frac{{\varphi }_1{\mu }_1{C}_{t1}}{k_1}{\frac{\partial m_1}{\partial t}}^{}\hspace{1em}{r}_w\le r<r_f$$

(3)

$$\frac{{\partial }^2{m}_2}{\partial r^2}+\frac{1}{r}\frac{\partial m_2}{\partial r}=\frac{{\varphi }_2{\mu }_2{C}_{t2}}{k_2}{\frac{\partial m_2}{\partial t}}^{}\hspace{1em}{r}_f\le r<\infty$$

(4)

where $r$ is the radial distance, $t$ is the time, $r_w$ is the wellbore radius, $k$ is the reservoir permeability, $\varphi$ is the porosity, $C_t$ is the total compressibility, and $r_f$ is the inner region radius (the subscripts 1 and 2 represent inner and outer regions, respectively). The initial condition for UGS is:

$$m_1\left(r,t=0\right)=m_2\left(r,t=0\right)=m_i$$

(5)

where $m_i$ is the initial pseudo−pressure. The outer boundary condition for UGS can be expressed as:

$$m_2\left(r\to \infty ,t\right)=m_i$$

(6)

The pressure and flow rate conditions for the interface between the inner and outer regions are:

$$m_1\left(r=r_f,t\right)=m_2\left(r=r_f,t\right)$$

(7)

$${\frac{\partial m_1}{\partial r}|}_{r=r_f}=\frac{k_2{\mu }_1}{k_1{\mu }_2}{\frac{\partial m_2}{\partial r}|}_{r=r_f}$$

(8)

The vertical and horizontal wells in UGS are considered to have constant rate production:

$${r\frac{\partial m_1}{\partial r}|}_{r=r_{\mathrm{w}}}=\frac{3.684\times {10}^{-3}{q}_{\mathrm{sc}}}{k_{1}h}\frac{p_{sc}T}{T_{sc}}$$

(9)

where $h$ is the reservoir thickness, $q_{sc}$ is the production rate in standard condition, $p_{sc}$ is the pressure in standard condition, and $T_{sc}$ is the temperature in standard condition. According to the dimensionless variables in Appendix A, the bottom−hole pressure solution of the vertical well can be expressed as Equation (10) (see Appendix B for details),

$${\overline{m}}_{wD}=a{I}_0\left(\sqrt{u}\right)+b{K}_0\left(\sqrt{u}\right)$$

(10)

where ${\overline{m}}_{wD}$ is the dimensionless bottom−hole pseudo−pressure, $u$ is the Laplace variable, $I_0$ is the zero−order first−class Bessel function, and $K_0$ is the zero−order second−class Bessel function. The undetermined coefficient $a$ and $b$ in Equation (10), be given as Equations (11) and (12):

$$a=\frac{m}{-u\sqrt{u}{I}_1(\sqrt{u})m+u\sqrt{u}{K}_1(\sqrt{u})n}$$

(11)

$$b=\frac{n}{-u\sqrt{u}{I}_1(\sqrt{u})m+u\sqrt{u}{K}_1(\sqrt{u})n}$$

(12)

Similarly, the undetermined coefficients $m$ and $n$ of Equations (11) and (12) can be expressed as:

$$m=K_0(r_{fD}\sqrt{u})\cdot \sqrt{\frac{u}{\sigma }}M{K}_1\left(r_{fD}\sqrt{\frac{u}{\sigma }}\right)-\sqrt{u}{K}_1(r_{fD}\sqrt{u})\cdot K_0\left(r_{fD}\sqrt{\frac{u}{\sigma }}\right)$$

(13)

$$n=-\sqrt{u}{I}_1(r_{fD}\sqrt{u})\cdot K_0\left(r_{fD}\sqrt{\frac{u}{\sigma }}\right)-I_0(r_{fD}\sqrt{u})\cdot \sqrt{\frac{u}{\sigma }}M{K}_1\left(r_{fD}\sqrt{\frac{u}{\sigma }}\right)$$

(14)

where $M$ is the diffusion ratio, $\sigma$ is the dispersion ratio, $r_{fD}$ is the dimensionless inner region radius, $I_1$ is the first−order first−class Bessel function, and $K_1$ is the first−order second−class Bessel function. The line source function method proposed by Ozkan was used to obtain the general pressure solution of horizontal wells in UGS [49]:

$${\overline{m}}_{wD}\left(x_D,y_D,u\right)=\frac{1}{2u}\left\{\begin{array}{c}{\displaystyle {\int }_{-1}^1\left[c_0{I}_0\left(\sqrt{u}{R}_D\right)+K_0\left(\sqrt{u}{R}_D\right)\right]d\varsigma }\\ +2{\displaystyle \sum _{n=1}^{\infty }{\displaystyle {\int }_{-1}^1\left[c_n{I}_0\left(a_n{R}_D\right)+K_0\left(a_n{R}_D\right)\right]d\varsigma }}\cos \left(n\pi z_D\right)\cos \left(n\pi z_{wD}\right)\end{array}\right\}$$

(15)

where $x_D$ and $y_D$ are the dimensionless horizontal distance, $z_D$ is the dimensionless vertical distance, $z_{wD}$ refers to the dimensionless vertical position for horizontal well, and $\varsigma$ is the integral variable. The $R_D$, $a_n$, $b_n$, and $c_n$ in Equation (15) can be given as Equations (16)–(19):

$$R_D=\sqrt{{\left(x_D-\varsigma \right)}^2+y_D{}^2}$$

(16)

$$a_n=\sqrt{u+n^2{\pi }^2{L}^2{}_{1D}}$$

(17)

$$b_n=\sqrt{\frac{u}{\sigma }+n^2{\pi }^2{L}^2{}_{2D}}$$

(18)

$$c_n=\frac{-M{b}_n{K}_0\left(a_n{r}_{fD}\right)K_1\left(b_n{r}_{fD}\right)+a_n{K}_1\left(a_n{r}_{fD}\right)K_0\left(b_n{r}_{fD}\right)}{M{b}_n{I}_0\left(a_n{r}_{fD}\right)K_1\left(b_n{r}_{fD}\right)+a_n{I}_1\left(a_n{r}_{fD}\right)K_0\left(b_n{r}_{fD}\right)}\begin{array}{cc}& n=0,1,2\dots \end{array}$$

(19)

The method proposed by Mukherjee and Economides was used to consider the effects of the wellbore storage coefficient and skin factor [44].

$${\overline{m}}_{wD}\left(C_D,S,u\right)=\frac{S+u{\overline{m}}_{wD}}{u\left[1+C_{D}u\left(u{\overline{m}}_{wD}+S\right)\right]}$$

(20)

where $C_D$ is the wellbore storage coefficient and $S$ is the skin factor. As given in Figure 3, the Stehfest numerical inversion method was used to obtain a real−domain bottom−hole pressure solution [45]. Figure 3 shows the solution procedure for vertical and horizontal wells in the gas storage.

2.3. Model Validation

In order to verify the correctness of the proposed model, the commercial numerical simulation software KAPPA was used to simulate the gas injection process for the vertical and horizontal wells in UGS, respectively. The reservoir domain included two regions, and the boundary for the inner region was circular. The outer boundary was chosen as infinite. The inner region radiuses for vertical and horizontal wells were chosen as 100 m and 500 m. The horizontal well’s length was 300 m. The mobility ratio and dispersion ratio between the inner and outer regions were set to 4. The wellbore storage factor and skin factor were selected as 0.5 m³/MPa and 2, respectively. In order to accurately describe the pressure change in the near−wellbore area, the minimum iteration time step was 1 × 10⁻⁵ h. The time step was gradually increased with the time increase to reduce the computational expense, and the maximum time step was 5 × 10⁴ h. The pressure and pressure derivatives from vertical and horizontal well models are given in Figure 4. There is a good agreement between the numerical results and the results from the proposed model. The input parameters for the model are given in

Table 1. Model validation parameters for vertical and horizontal wells with composite features in underground natural gas storage.

Parameter	Value	Unit
UGS thickness	25	m
Permeability	35	mD
porosity	0.21	−
Initial pressure	34.478	MPa
Total compressibility	0.0192	MPa−1
Production rate	80 × 104	m3/d
Wellbore storage coefficient	0.5	m3/MPa
Skin factor	2	−
Horizontal well length	300	m
Composite distance of the vertical well	100	m
Composite distance of the horizontal well	500	m
Mobility ratio	4	−
Dispersion ratio	4	−

.

3. Flow Chart

With the theoretical method in Section 2.2, the storage performance evaluation technical process is shown in Figure 5. Specifically, the following main steps can be included:

1.

Step 1 Inventory Parameters Evaluated by Gas Reservoir Engineering

Lessons from the development of gas reservoirs indicate that the inventory evaluation parameters include inventory, effective inventory, effective gas volume, peak shaving gas volume, working gas volume, and cushion gas volume [50,51]. Similar to SEC reserves, effective inventory, effective gas volume, and working gas volume are key evaluation parameters for listing transactions, and they directly determine the market valuation of UGS companies. As shown in Figure 6, the effective inventory is equivalent to the technically recoverable reserves in gas reservoir engineering, and it is one of the parameters to measure the UGS performance. The material balance method in Equation (21) proposed by Tek can be used to calculate the effective inventory value [41].

$$G=\frac{Q_t}{{\left(\frac{p}{ZT}\right)}_{gi}-{\left(\frac{p}{ZT}\right)}_{gp}}{\left(\frac{p}{ZT}\right)}_{gi}$$

(21)

where $G$ is the effective inventory, $Q_t$ is the total amount of produced gas, the subscript $gi$ is the ending time of gas injection stage, and $gp$ is the ending time of gas withdrawal stage. The effective gas volume directly determines the gas storage capacity of UGS. According to the gas state equation, the effective gas volume can be calculated by Equation (22).

$$V=G\frac{Z_{gi}{T}_{gi}}{p_{{}_{gi}}}\frac{p_{sc}}{T_{sc}}$$

(22)

where $V$ is the effective gas storage volume, $p_{{}_{gi}}$ is the AFP at the end of the gas injection stage, $Z_{{}_{gi}}$ is the deviation factor at the end of the gas injection stage, and $T_{{}_{gi}}$ is the formation temperature. The working gas volume is used to meet market demand. As shown in Figure 6, the working gas volume can be calculated using Equation (23).

$$G_w=V\frac{p_{\max }}{Z_{\max }{T}_{gi}}\frac{T_{sc}}{p_{sc}}-V\frac{p_{\min }}{Z_{\min }{T}_{gi}}\frac{T_{sc}}{p_{sc}}$$

(23)

where $G_w$ is the working gas volume, $p_{\max }$ is the upper limit pressure, $p_{\min }$ is the lower limit pressure, and $Z_{\max }$ and $Z_{\min }$ are deviation factors for the upper limit pressure and the lower limit pressure, respectfully.

2.

Step 2: Error Sources in Error Theory

Due to the limitations of monitoring equipment and calculation methods, there will be errors in the monitoring results of temperature, pressure, gas injection, and withdrawal volume in UGS. With the error theory proposed by Tek, as shown in Equation (24), the error sources in UGS mainly come from the AFP, deviation factor, pseudo−pressure, and pseudo−pressure drop [41]. In addition, deviation factor, pseudo−pressure, and pseudo−pressure drop are related to the formation pressure. Thus, the AFP plays a more critical role in inventory evaluation. The detailed derivation for Equation (24) can be found in Appendix C.

$$\epsilon (G)=\frac{\mathrm{\Delta }{p}_{gi}}{p_{gi}}-\frac{\mathrm{\Delta }{Z}_{gi}}{Z_{gi}}-\frac{{\left(\frac{p}{ZT}\right)}_{gi}}{{\left(\frac{p}{ZT}\right)}_{gi}-{\left(\frac{p}{ZT}\right)}_{gp}}\left(\frac{\mathrm{\Delta }{p}_{gi}}{p_{gi}}-\frac{\mathrm{\Delta }{Z}_{gi}}{Z_{gi}}\right)-\frac{{\left(\frac{p}{ZT}\right)}_{gp}}{{\left(\frac{p}{ZT}\right)}_{gi}-{\left(\frac{p}{ZT}\right)}_{gp}}\left(\frac{\mathrm{\Delta }{p}_{gp}}{p_{gp}}-\frac{\mathrm{\Delta }{Z}_{gp}}{Z_{gp}}\right)$$

(24)

where $\epsilon$ is the relative error and $\mathrm{\Delta }$ is the absolute error.

3.

Step 3: Calculation of Average Formation Pressure

Different from the gas reservoir development process, the gas injection and withdrawal volume in UGS is a constant value that is determined by social demand. Due to the requirements of society’s heating needs, the gas injection phase in UGS usually occurs from April to November, and the gas withdrawal process usually includes January, February, and March. This long−term and large−scale gas injection and withdrawal process leads to uneven pressure distribution in UGS, and the well interference is serious [52,53,54]. We use the radial composite model and the homogeneous model to calculate the AFP for the single well. However, the pressure for a single well often cannot represent the AFP of the entire UGS. The pressure superposition principle and weighting method of the gas injection−withdrawal rate are used to calculate the AFP. Since the target partial differential equation is linear, the pressure superposition principle in Figure 7 can be used to obtain the pressure solution with the effect of adjacent wells through Equations (10), (15) and (20) [55,56]. For wells that are farther away, the weighting method of gas injection−withdrawal rate is used [42].

$${\overline{m}}_{wD}(u)={\overline{m}}_{iwD}(u)+{\displaystyle \sum _{j=2}^{N_w}{\overline{m}}_{j,D}(r_{j,D},u)},j=2,3\dots N_w$$

(25)

where ${\overline{m}}_{iwD}$ is the bottom−hole pressure solution of the target well, ${\overline{m}}_{j,D}$ is the pressure interference from the adjacent well, $r_{j,D}$ is the dimensionless distance between the adjacent well and the target well, and $N_w$ is the number of wells in UGS. The weighting method of the gas injection−withdrawal rate can be described by Equation (26).

$$p_a={\displaystyle \sum _{k=1}^{N_w}{p}_k\frac{Q_k}{{\displaystyle \sum _{j=1}^N{Q}_j}}}$$

(26)

where $p_a$ is the AFP of gas storage, $Q$ is the cumulative gas injection and production volume of a single well, and the subscripts $k$ and $j$ represent serial numbers of wells.

4.

Step 4: Calculation of the Deviation Factor

The deviation factor reflects the deviation degree between the real gas and the ideal gas, and it is used to correct the ideal gas state equation [57]. The method proposed by the American Gas Association (AGA) is used to calculate the deviation factor [58]. This method requires a detailed molar composition analysis of the gas mixture. The components with a mole fraction less than 0.00005 are ignored. In addition, hydrogen sulfide, water, and hydrocarbon components with carbon molecular number greater than 8 cannot be considered. It can be expressed as Equation (27) (please see Appendix D for details).

$$Z=1+B{\rho }_m-{\rho }_r{\displaystyle \sum _{n=13}^{18}{C}_n^{*}}+{\displaystyle \sum _{n=13}^{58}{C}_n^{*}}(b_n-c_n{k}_n{\rho }_r^{k_n}){\rho }_r^{b_n}\exp (-c_n{\rho }_r^{k_n})$$

(27)

where ${\rho }_r$ is the reduced density, $b_n$, $c_n$, and $k_n$ are undetermined coefficients, $B$ is the second virial coefficient, and $C_n^{*}$ is the equation coefficient related to gas composition and temperature. Combining Equations (2) and (27), we can get the final equation:

$$p={\rho }_{m}RT\left[1+B{\rho }_m-{\rho }_r{\displaystyle \sum _{n=13}^{18}{C}_n^{*}}+{\displaystyle \sum _{n=13}^{58}{C}_n^{*}}(b_n-c_n{k}_n{\rho }_r^{k_n}){\rho }_r^{b_n}\exp (-c_n{\rho }_r^{k_n})\right]$$

(28)

Equation (28) is the transcendental equation in one variable about ${\rho }_m$. The bisection method is used to calculate the pressure value. When the absolute error of pressure is less than 1 × 10⁻⁶, the calculation process is considered to be finished.

4. Field Application

This section uses the actual data of pressure and gas injections, as well as the production volumes of Hutubi UGS to carry out a field application. First, the geological background of Hutubi UGS and its reservoir properties, fluid properties, and well properties are briefly introduced. The typical pressure monitoring results for Hutubi UGS are discussed. The storage inventory parameters from Hutubi UGS are determined using the proposed method.

4.1. Geological Background

The Hutubi UGS is China’s first UGS for receiving gas imported from Central Asia, and it is located at the eastern end of the third row of structural belts in the North Tianshan piedmont depression of the Junggar Basin. The gas reservoir is a set of retrogradational delta deposits in a lacustrine setting, which is a large continuous sedimentary depression [59]. The gas storage layer is the Tertiary Ziniquanzi formation, including the Zi I and Zi II layers. The top depth of the UGS is about 3500 m. The thicknesses of the gas reservoir, direct caprock, and regional caprock are 120 m, 150 m, and 870 m, respectively [60]. The reservoir lithology is mainly fine sandstone and siltstone, and the pores are primary intergranular pores. The average porosities for Zi I and Zi II layers are 19.4% and 15.4%, respectively. The absolute permeabilities for two layers are 48.6 mD and 19.7 mD. The Hutubi gas field started to produce in April 1998. In 2010, the gas production from the Hutubi gas reservoir was below the economic production value, and it started the construction of the UGS. The remaining geological reserve before the UGS construction was 45.3 × 10⁸ m³ [61]. The UGS has been in operation since June 2013. Currently, the Hutubi UGS includes 29 vertical wells and 14 horizontal wells. Gas chromatographic data shows that the methane content in the produced gas from the Hutubi UGS is above 92.5% and the content of non−hydrocarbon components is less than 2%. The average relative density of injected gas is 0.59. The designed upper limit pressure and lower limit pressure for the Hutubi UGS are 34 MPa and 18 MPa, respectively. The designed storage capacity, working gas volume, and maximum daily gas injection volume are 107 × 10⁸ m³, 45.1 × 10⁸ m³, and 1550 × 10⁴ m³ [62,63,64]. The geological structure and well location map of the Hutubi UGS are shown in Figure 8.

4.2. Pressure Monitoring

Three pressure monitoring methods have been applied to the Hutubi UGS, namely the permanent monitoring well, the direct monitoring for bottom−hole pressure, and the pressure transient monitoring [61,67]. The permanent monitoring wells in the Hutubi UGS include one well for caprock monitoring, nine wells for trap, and three wells for formation. The trap monitoring results show that the pressure fluctuates between 28.7 and 31.3 MPa. The small variation in pressure indicates that the gas storage is in a sealed state. The monitoring pressure inside the formation varies widely and is injection−withdrawal cycle−related. The pressure range during the 1–4th injection−withdrawal cycle is 14.6–30.8 MPa. Subsequently, the pressure tends to be stable, and the range is 22.3–32.9 MPa. In addition, the bottom−hole pressure of the injection−withdrawal well after 15 days of the shut−in period is selected as the direct monitoring result for bottom−hole pressure. The trend of pressure change is similar to the results of permanently monitored wells. The pressure range for the first four gas injection cycles is 19.9–34.0 MPa. For gas production duration, the range is 18.4–26.6 MPa. After the fourth injection−withdrawal cycle, the above two pressure ranges were updated to 29.5–34.2 MPa and 23.0–26.7 MPa, respectively. To evaluate the injection−production capacity of gas storage wells and the AFP, 127 pressure transient tests are also carried out during nine injection−withdrawal cycles. Gas injection and withdrawal testing times are 210 h and 240 h, respectively. The obtained testing pressure data were matched by the proposed model, and the balanced pressure values were obtained after a prolonged testing time. Compared with the direct monitoring results, the balanced pressure ranges during the 1–4th injection−withdrawal cycle changed by 4.3% and 3.5%, respectively. After the fourth injection−withdrawal cycle, the pressure ranges for gas injection and withdrawal are 2.3% and 2.2%, respectively. The transient pressure behaviors in Figure 9 show that the potential flow regimes in the Hutubi UGS include wellbore storage effect flow, skin effect flow, linear flow (for horizontal wells), pseudo−radial flow in the inner zone, transitional flow, and pseudo−radial flow in the outer zone. Different from the development of gas reservoirs, the large−scale injection−withdrawal rate in UGS leads to the obvious radial composite characteristics in pressure data. With the increase of the injection−withdrawal cycle, the derivative curve gradually rises, and the radial composite feature becomes obvious. The evaluated permeability distribution range is 16–55 mD. The skin factor and mobility ratio ranges are 0–5 and 2–6.2. The inner zone radius distribution ranges from 50 m to 300 m.

4.3. Model Application

With the proposed methodology, the inventory evaluation for the Hutubi UGS is determined by the following steps.

1.

Step 1: Inventory Parameters Evaluated by Gas Reservoir Engineering

Based on the gas reservoir engineering method, the material balance equation in Equation (21) and the gas state equation in Equations (22) and (23), the main technical parameters for UGS inventory evaluation include inventory, effective inventory, effective gas volume, and working gas volume.

2.

Step 2: Error Sources in Error Theory

The relative error of effective inventory in Equation (24) shows that the error sources in inventory evaluation are mainly the relative errors of AFP, deviation factor, pseudo−pressure, and pseudo−pressure drop. Deviation factor, pseudo−pressure, and pseudo−pressure drop are also related to the AFP. The AFP is therefore the main source of inventory errors.

3.

Step 3: Calculation of the AFP and Deviation Factor

First, we collected the reservoir, fluid, and well properties of the Hutubi UGS in order to calculate the AFP. Subsequently, we collected the transient pressure testing data of 127 wells, the direct pressure monitoring data of 415 wells, and the permanent monitoring data of 13 wells from the Hutubi UGS. With curve matching of pressure transient testing data, the typical parameters, including permeability, skin coefficient, mobility ratio, and inner zone radius, can be obtained. The pressure solutions from vertical and horizontal wells in Equations (10) and (15) are obtained by the Laplace transform technology. The wellbore storage and skin effects are considered in the method of Mukherjee and Economides [44]. The balanced state pressure of a single well in the gas storage is obtained by prolonging the well shut−in time. The pressure superposition method in Equation (25) and the weighting method of the gas injection−withdrawal rate in Equation (26) are used to obtain the AFP. The weighting coefficient is related to the gas injection and production volume, and the weighting coefficient of Hutubi UGS is between 0.02 and 0.12. As shown in Figure 10, the balanced AFP for the Zi I layer at the end of the injection−withdrawal cycle increased rapidly from 18.7 MPa and 16.3 MPa in the 1st cycle to 32.1 MPa and 25.4 MPa in the 4th cycle. After the fourth period, the AFP tends to remain unchanged. The Zi I layer AFP for the injection−withdrawal cycle is finally stable at around 32.9 MPa and 22.2 MPa. For the Zi II layer, the final AFP values are 32.8 MPa and 22.4 MPa. The reason is that the gas injection volume of UGS at the initial stage is greater than the production volume to replenish formation energy. According to the AFP, for reservoir temperature and gas composition,

Table 2. The deviation factor calculation results for Hutubi UGS during nine injection−withdrawal cycles.

Injection−Withdrawal Cycle	Gas Injection Cycle		Gas Withdrawal Cycle
	Zi I Layer	Zi II Layer	Zi I Layer	Zi II Layer
1	0.933	0.930	0.930	0.929
2	0.956	0.932	0.941	0.931
3	0.989	0.944	0.953	0.937
4	1.007	0.972	0.960	0.952
5	1.006	0.994	0.960	0.955
6	1.006	1.001	0.964	0.960
7	1.002	1.002	0.963	0.959
8	1.019	1.016	0.948	0.950
9	1.014	1.013	0.944	0.945

shows that the deviation factor in Equations (2), (27) and (28) varies from 0.93 to 1.02.

4.

Step 4: Inventory Evaluation Results

After completing the above steps, the evaluation parameters, including effective inventory, effective gas volume, and working gas volume, are obtained by Equations (21)–(23). The effective inventory, effective storage volume, and working gas volume in Figure 11 increase rapidly at the initial stage. After entering the fourth injection−production cycle, the increase rate of effective inventory, effective gas volume, and working gas volume slows down. After nine injection−withdrawal cycles, Figure 12 shows that the inventory of Hutubi UGS has increased to 100.1 × 10⁸ m³ and the capacity filling rate is 93.6%. The inventory utilization ratio increased from 41.9% to 95.2%. The final effective inventory, effective gas volume, and working gas volume are 95.3 × 10⁸ m³, 3626 × 10⁴ m³, and 40.3 × 10⁸ m³.

5. Summary and Conclusions

This study developed a novel method for performance evaluation of the gas reservoir−type UGS. Subsequently, the proposed workflow was applied to the field case of the largest Hutubi UGS in China. Some obtained key conclusions are as follows:

• The workflow and theoretical basis behind the method were first introduced. The flexibility and effectiveness of this method in UGS performance evaluation were illustrated. The well models based on the Laplace transform technology and the Stehfest numerical inversion method were used to match the field pressure monitoring data. By prolonging the well shut−in time, the balanced formation pressure represented by a single−well was obtained. The pressure superposition principle and weighting method of the gas injection−withdrawal rate were selected to obtain the AFP.
• The pressure monitoring results show that large−scale injection and production lead to a continuous rise of the pressure derivative curve in the middle and late stages. This reflects the formation heterogeneity caused by gas injection and production. For the Zi I layer, the AFP at the end of the gas injection and withdrawal cycle increased rapidly from 18.7 MPa and 16.3 MPa in the first cycle to 32.1 MPa and 25.4 MPa in the fourth cycle. After the fourth cycle, the AFP tended to be stable. The final AFP values of the Zi I layer were stable at around 32.9 MPa and 22.2 MPa. For the Zi II layer, the final AFP values were 32.8 MPa and 22.4 MPa.
• By using the material balance method, the inventory, effective inventory, and working gas volume were calculated. During the nine injection−withdrawal cycles in Hutubi UGS, the changes in inventory parameters could be divided into a rapid−increase stage and a steady stage. In the rapid−increase stage, the inventory, effective inventory, and working gas volume increased to 97.7 × 10⁸ m³, 84.5 × 10⁸ m³, and 35.7 × 10⁸ m³, respectively, in the fourth cycle. From the fourth cycle, these three parameters began to stabilize, and the variation range was less than 11.4%. The final inventory, effective inventory, and working gas volume were 100.1 × 10⁸ m³, 95.3 × 10⁸ m³, and 40.3 × 10⁸ m³, respectively.