Utilization of CO 2 as Cushion Gas for Depleted Gas Reservoir Transformed Gas Storage Reservoir

: Underground gas storage reservoirs (UGSRs) are used to keep the natural gas supply smooth. Native natural gas is commonly used as cushion gas to maintain the reservoir pressure and cannot be extracted in the depleted gas reservoir transformed UGSR, which leads to wasting huge amounts of this natural energy resource. CO 2 is an alternative gas to avoid this particular issue. However, the mixing of CO 2 and CH 4 in the UGSR challenges the application of CO 2 as cushion gas. In this work, the Donghae gas reservoir is used to investigate the suitability of using CO 2 as cushion gas in depleted gas reservoir transformed UGSR. The impact of the geological and engineering parameters, including the CO 2 fraction for cushion gas, reservoir temperature, reservoir permeability, residual water and production rate, on the reservoir pressure, gas mixing behavior, and CO 2 production are analyzed detailly based on the 15 years cyclic gas injection and production. The results showed that the maximum accepted CO 2 concentration for cushion gas is 9% under the condition of production and injection for 120 d and 180 d in a production cycle at a rate of 4.05 kg / s and 2.7 kg / s, respectively. The typical curve of the mixing zone thickness can be divided into four stages, which include the increasing stage, the smooth stage, the suddenly increasing stage, and the periodic change stage. In the periodic change stage, the mixed zone increases with the increasing of CO 2 fraction, temperature, production rate, and the decreasing of permeability and water saturation. The CO 2 fraction in cushion gas, reservoir permeability, and production rate have a signiﬁcant e ﬀ ect on the breakthrough of CO 2 in the production well, while the e ﬀ ect of water saturation and temperature is limited.


Introduction
Natural gas storage is principally used for meeting the widely fluctuating demand of gas, especially the high peak demands in winter [1][2][3]. Meanwhile, the development of natural gas storage is further promoted by a profitable business model, i.e., storing at a lower price and selling at a higher price based on demand loads. [4]. There are several vital strategies for natural gas storage such as gas tanks, salt caverns, and gas pipeline [5]. Owning to its advantages of a large storage capacity and storage is mostly a seasonal operation that usually is charged in summer and discharged in winter. In the operation process, cushion gas is vital to maintain suitable reservoir pressure as well as to keep a stable production operation.
In the case of depleted gas reservoir transformed gas storage reservoir, the native natural gas is commonly used as the cushioning gas [12]. This method may lead to a great deal of waste because approximately 40%-70% of stored gas works as cushion gas to provide pressure support and cannot be extracted [13]. In addition, the security of a UGSR during the CH4 injection could be threatened dramatically by the increased reservoir pressure due to the relatively low compressibility of the cushion gas. The utilization of supercritical fluids as cushion gas is a potential strategy to tackle these issues. For supercritical fluids there is a Widom region, where some physicochemical properties such as density and compressibility exhibit anomalous behavior [14][15][16]. The Widom region is related to the critical pressure and critical temperature, which is 7.38 MPa and 31.04 °C, respectively for CO2. This critical condition exactly can be achieved in the storage applications, demonstrating the potential of utilizing CO2 as cushion gas in the natural gas storage reservoir. In comparison with the native gas cushion, 30% more CH4 can be stored with CO2 as cushion gas [12]. It may also bring an additional benefit to enhance gas recovery in the post-CO2 storage and enhanced gas recovery (CSEGR) process for a depleted gas reservoir. However, the mixing of CO2 and CH4 with the injection and production of working gas as shown in Figure 1 can be an obstacle for the utilization of CO2 as cushion gas [12], which requires serious attention. Some research efforts have been paid to study the mixing behavior of CO2 and CH4 in UGSR. Oldenburg [12] conducted a numerical simulation on a two-dimensional reservoir model, to demonstrate the suitability of CO2 as cushion gas in natural gas storage. Based on this reservoir model, Ma et al. [17] studied the impact of geological parameters on the gas mixing behavior in UGSR with CO2-based cushion gas by using a hydromechanically coupled model. Their results showed that the mixing region decreases with increasing reservoir thickness and dip angle. However, the reservoir model has been initially saturated with CO2, which neglects the replacing process of the native CH4 by the CO2 during its injection into the depleted gas reservoir. In addition, it did not characterize the seasonal injection and production, as the mixing behavior has been studied only within 180 days. Moreover, the presence of residual water in the reservoir, which could affect the  Some research efforts have been paid to study the mixing behavior of CO 2 and CH 4 in UGSR. Oldenburg [12] conducted a numerical simulation on a two-dimensional reservoir model, to demonstrate the suitability of CO 2 as cushion gas in natural gas storage. Based on this reservoir model, Ma et al. [17] studied the impact of geological parameters on the gas mixing behavior in UGSR with CO 2 -based cushion gas by using a hydromechanically coupled model. Their results showed that the mixing region decreases with increasing reservoir thickness and dip angle. However, the Energies 2020, 13,576 3 of 20 reservoir model has been initially saturated with CO 2 , which neglects the replacing process of the native CH 4 by the CO 2 during its injection into the depleted gas reservoir. In addition, it did not characterize the seasonal injection and production, as the mixing behavior has been studied only within 180 days. Moreover, the presence of residual water in the reservoir, which could affect the mixing behavior, has not been considered in their work. Niu and Tan [18] investigated the impacts of reservoir porosity and initial operation pressure on the mixing behavior of the cushion and working gases based on the three-dimensional gas-water two phase theory. Their results showed that the reservoir with high porosity and a large initial pressure is favored for limiting the formation and migration of the mixed zone. Oldenburg and Pan [19] conducted a one-dimensional radial simulation to investigate the pressurization and gas-gas mixing behavior. It was found that the pressure rise during working gas injection reduces if CO 2 is used as cushion gas. They also discovered that the impact of the CO 2 cushion on the reservoir pressure during the production period is much more significant than the one during the injection period, due to a higher production rate compared to the injection rate. The behavior of CO 2 and N 2 as cushion gases were compared by Kim et al. [20], the results showed that CO 2 is more suitable as cushion gas in depleted gas reservoirs in terms of the productivity index. However, in these works, a single well for both cushion gas injection and working gas production is very likely to accentuate the mixing of gases in a UGSR, which is not favorable in the engineering field.
For the purpose of investigating the suitability of using CO 2 as the cushion gas, the Donghae depleted gas reservoir located in Ulleung basin in Korea is used. It is important to note that the residual water is also considered in this work. After the injection of CO 2 , the cyclic CH 4 production and injection was conducted over a period of 15 years, which is closer to a real engineering operation and makes it more beneficial for understanding the mixing behavior of CO 2 and CH 4 in a relatively long-term period. Additionally, the effect that the CO 2 fraction, reservoir temperature, reservoir permeability, residual water, and production rate have on the mixing behavior and gas production were analyzed in detail.

Governing Equations
The multicomponent and multiphase flow TOUGH2MP/TMVOC simulator, which could characterize the behavior of CO 2 and CH 4 in unsaturated zones, was used in this study [21]. The basic mass and energy balance equations for multi-component and multi-phase fluid flow in porous medium can be written in following form [22]: where V n is the control volume for an arbitrary subdomain in the flow system; Γ n represents the closed surface of V n ; n is a normal vector pointing inward into V n ; M k denotes the mass or energy per volume; k is the mass or energy components; F k denotes the mass or heat flux; and q k represents the sinks and sources [22]. The mass accumulation term for water, CO 2 and CH 4 in the systems can be expressed as: where φ is the porosity; S β and ρ β are the saturation and density of phase β, respectively; and X k β is the mass fraction of component k in phase β. The advective mass flux F k is a sum of phases which is given by: Energies 2020, 13, 576 4 of 20 The individual phase F β is calculated by the multiphase version of Darcy's law: where u β is the Darcy velocity in phase β; K is the absolute permeability; K rβ denotes the relative permeability to phase β; µ β is viscosity; g denotes the vector of gravitational acceleration; and P β represents the fluid pressure in phase β, which can be written as: where P is the pressure of reference phase and P cβ is the capillary pressure. The diffusive flux of water, CO 2 and CH 4 in phase β are given in: where τ 0 is a porous medium dependent factor; τ β is a coefficient determined by phase saturation; d k β denotes the diffusion coefficient of component k in phase β; and x k β represents the mole fraction of component k in phase β.
A single effective multiphase diffusion coefficient can be defined as: For two-phase conditions in the UGSR, total diffusive flux is then expressed by: The pressure and temperature dependent diffusion coefficients for CO 2 and CH 4 is given by: where P 0 and T 0 denote the standard conditions, which equal 0.1 MPa and 0 • C, respectively; the temperature dependence parameter θ is 1.8. The error percentage of the diffusion coefficients for CO 2 and CH 4 calculated form Equation (9) is between 2.31% and 12.78% compared with the experimental results of Honari et al. [23] for the Estaillades carbonate, Donnybrook sandstone, and Ketton carbonate, demonstrating the suitability of this empirical model. This model was implemented into the multicomponent and multiphase flow simulator TOUGH2MP that was used in this work.

Geological Model
A simple brick model with the dimension of 914.4 m × 914.4 m × 30.48 m was generated based on the geological information of the Donghae gas reservoir from the Ulleung basin in Korea, in which a slightly anticline structure is neglected [20]. As shown in Figure 2a, a typical five-spot pattern was used in this study, wherein one well is located at the lower formation for CO 2 injection, and the other four wells are located at upper formation for the CH 4 injection and production. In this symmetrical model, only a quarter model of the reservoir is selected for the simulation (Figure 2b). The applied geological properties of the reservoir are summarized in Table 1. on the geological information of the Donghae gas reservoir from the Ulleung basin in Korea, in which a slightly anticline structure is neglected [20]. As shown in Figure 2a, a typical five-spot pattern was used in this study, wherein one well is located at the lower formation for CO2 injection, and the other four wells are located at upper formation for the CH4 injection and production. In this symmetrical model, only a quarter model of the reservoir is selected for the simulation (Figure 1b). The applied geological properties of the reservoir are summarized in Table 1.

Operation Scenarios
Considering the low pressure of the depleted reservoir, two years of gas injection, i.e., CO 2 injection in the first year and CH 4 injection in the second year, was implemented before its operation as a gas storage reservoir, to recover the average reservoir pressure from 5.17 MPa to the original reservoir pressure of 24 MPa [24]. As shown in Figure 3, the gas storage reservoir works in one-year cycles and each cycle consists of four stages. In the first stage, CH 4 was produced at a rate of 4.05 kg/s from 1 November to 28 February the next year (120 days), followed by well shutting and facility checking from 1 March to 4 April (35 days). In the third stage, CH 4 was injected at a rate of 2.7 kg/s from 5 April to 1 October (180 days). Then the well was shut again from 2 October to 31 October (30 days).

Operation Scenarios
Considering the low pressure of the depleted reservoir, two years of gas injection, i.e., CO2 injection in the first year and CH4 injection in the second year, was implemented before its operation as a gas storage reservoir, to recover the average reservoir pressure from 5.17 MPa to the original reservoir pressure of 24 MPa [24]. As shown in Figure 3, the gas storage reservoir works in one-year cycles and each cycle consists of four stages. In the first stage, CH4 was produced at a rate of 4.05 kg/s from 1 November to 28 February the next year (120 days), followed by well shutting and facility checking from 1 March to 4 April (35 days). In the third stage, CH4 was injected at a rate of 2.7 kg/s from 5 April to 1 October (180 days). Then the well was shut again from 2 October to 31 October (30 days).

Physical Properties of the Mixed Gases
The mixing of CO2 and CH4 in a UGSR is affected by many factors, such as density differences, mobility ratios, molecular diffusion, and mechanical dispersion [2]. The density difference between

Physical Properties of the Mixed Gases
The mixing of CO 2 and CH 4 in a UGSR is affected by many factors, such as density differences, mobility ratios, molecular diffusion, and mechanical dispersion [2]. The density difference between CO 2 and CH 4 plays the most important role in the separation of the gases. Figure 4a,b shows the density and viscosity of CO 2 -CH 4 mixtures at 40 • C, respectively, which were calculated by the WebGasEOS v.2.01 developed by the Lawrence Berkeley National Laboratory [25]. In comparison with CH 4 , the higher density of CO 2 could lead to downward sinking of CO 2 . As shown in Figure 4a, the density of the mixed gas was strongly correlated to the gas composition and pressure. Especially, with addition of small amounts the CH 4 into CO 2 , the density decreased sharply. Such decreasing behavior slowed down with increasing CH 4 fraction. It should be mentioned that the sudden increase of density occurred near the critical pressure of CO 2 because of the CO 2 entering into the supercritical state.
Energies 2020, 13, 576 6 of 21 behavior slowed down with increasing CH4 fraction. It should be mentioned that the sudden increase of density occurred near the critical pressure of CO2 because of the CO2 entering into the supercritical state. The mobility differences in CO2-CH4 displacement was primarily caused by different dynamic viscosities. Figure 4b shows the dynamic viscosities of CO2-CH4 mixtures, whose tendency is similar as that of the density in Figure 4a. The difference of viscosities between CO2 and CH4 is also favorable for limiting the gas mixing [26], but it may lead to an unstable contact interface [17].  The dispersion coefficient of CO2 in CH4 normally ranges from 0.01 to 0.3 cm 2 /min [27], which is a relatively slow velocity compared with the advective and convective transport. It should be pointed that the diffusion effect is proportional to concentration gradient, thus its effect will decrease along with the mixing process of gases. The mechanical dispersion is controlled by the movement of formation fluids [2].
3.1.2. Spatial Distribution of CO2 and CH4 Figure 5 shows the spatial distribution of CO2 at the end of the CH4 production stage, well shutting stage, CH4 injection stage, and the second well shutting stage, i.e., the time of 120 d, 155 d, 335 d, and 365 d in the 1st, 5th, 10th, and 15th year, respectively for the gas storage reservoir with  The mobility differences in CO 2 -CH 4 displacement was primarily caused by different dynamic viscosities. Figure 4b shows the dynamic viscosities of CO 2 -CH 4 mixtures, whose tendency is similar as that of the density in Figure 4a. The difference of viscosities between CO 2 and CH 4 is also favorable for limiting the gas mixing [26], but it may lead to an unstable contact interface [17]. The dispersion coefficient of CO 2 in CH 4 normally ranges from 0.01 to 0.3 cm 2 /min [27], which is a relatively slow velocity compared with the advective and convective transport. It should be pointed that the diffusion effect is proportional to concentration gradient, thus its effect will decrease along with the mixing process of gases. The mechanical dispersion is controlled by the movement of formation fluids [2].
3.1.2. Spatial Distribution of CO 2 and CH 4 Figure 5 shows the spatial distribution of CO 2 at the end of the CH 4 production stage, well shutting stage, CH 4 injection stage, and the second well shutting stage, i.e., the time of 120 d, 155 d, 335 d, and 365 d in the 1st, 5th, 10th, and 15th year, respectively for the gas storage reservoir with 10% CO 2 as the cushion gas. It can be seen that the CO 2 was concentrated at the bottom of the reservoir due to the large density compared with CH 4 . The concentration of CO 2 in the reservoir decreased gradually with time, resulting from the production of CO 2 and the mixing of CO 2 with CH 4 . Specifically, the maximum concentration of CO 2 in the UGSR decreased from 100% to 92.3%, 68.7%, 53.5%, and 43.6% at the end of the 1st, 5th, 10th, and 15th year, respectively. The distribution of CO2 changes periodically over time corresponding to the operation scenarios ( Figure 3). In the first stage (CH4 production), CO2 moved towards the production well due to the pressure gradient. Meanwhile, the gases in the reservoir were significantly mixed until the end of gas production. In the second stage (well shutting), the behavior of the mixed gases was controlled dominantly by the density difference and diffusion effect. However, the effect on mixing was still minor due to the relatively low diffusion and short well shut-in time. In the third stage (CH4 injection), CO2 migrated towards the CO2 injection well under the displacement effect of CH4 injection. The mixed zone of gases became smaller in this stage. The last stage, i.e., the second well shutting stage, also had a negligible effect on the gases mixing.  The distribution of CO 2 changes periodically over time corresponding to the operation scenarios ( Figure 3). In the first stage (CH 4 production), CO 2 moved towards the production well due to the pressure gradient. Meanwhile, the gases in the reservoir were significantly mixed until the end of gas production. In the second stage (well shutting), the behavior of the mixed gases was controlled dominantly by the density difference and diffusion effect. However, the effect on mixing was still minor due to the relatively low diffusion and short well shut-in time. In the third stage (CH 4 injection), CO 2 migrated towards the CO 2 injection well under the displacement effect of CH 4 injection. The mixed zone of gases became smaller in this stage. The last stage, i.e., the second well shutting stage, also had a negligible effect on the gases mixing.

Effect of the CO 2 Fraction
To figure out the suitability of CO 2 as cushion gas, the UGSR with a CO 2 concentration varying from 0% to 20% for the cushion was investigated. The reservoir average pressure against time is shown in Figure 6. It can be seen that all of the reservoir average pressure changed periodically over time, decreasing in the production stage, remaining constant during well shutting stage, and increasing again in the injection stage. In the first cycle, the minimum value of reservoir pressure decreased with the increasing CO 2 concentration due to the higher compressibility of CO 2 compared with CH 4 . With the operation of the UGSR, the maximum reservoir average pressure increased gradually as the produced CO 2 was replaced by the injected CH 4 , which had a lower compressibility, and the increase in CO 2 concentration led to pressure increment.

Effect of the CO2 Fraction
To figure out the suitability of CO2 as cushion gas, the UGSR with a CO2 concentration varying from 0% to 20% for the cushion was investigated. The reservoir average pressure against time is shown in Figure 6. It can be seen that all of the reservoir average pressure changed periodically over time, decreasing in the production stage, remaining constant during well shutting stage, and increasing again in the injection stage. In the first cycle, the minimum value of reservoir pressure decreased with the increasing CO2 concentration due to the higher compressibility of CO2 compared with CH4. With the operation of the UGSR, the maximum reservoir average pressure increased gradually as the produced CO2 was replaced by the injected CH4, which had a lower compressibility, and the increase in CO2 concentration led to pressure increment. To quantify the spatial extent of the mixing region, the mixed thickness, defined as the distances between the CO2 fraction level of 10% and 90% along the diagonal line crossing both the CO2 injection well and the CH4 production well (see in Figure 7), were compared in Figure 8 under various initial concentrations of CO2 the cushion gas. The mixed thickness was used to characterize the interface between the working and cushion gases [19]. Clearly, the thickness of the mixed zone increased with the increasing of the CO2 initiated concentration. The curve of the mixed zone's thickness for the UGSR with initial CO2 concentration ranging from 8% to 10% was similar, in which the representative curve for the initial 9% of CO2 could be divided into four stages over the lifetime of project (see in Figure 9c). Besides, the distance between the CO2 injection well and 10% CO2 (r0.1), a 90% CO2 (r0.9) concentration point along the same diagonal line is plotted in Figure 9a,b against time, respectively. The variation of r0.1 is harmonized with the periodical operation scenarios (Figure 3). It shows that r0.1 was mainly driven by the gas displacement effect. However, this value was limited in a narrow range from 189.9 to 193.7 m. Therefore, it could be concluded that the abrupt change of mixed zone thickness (r0.1 − r0.9) was mostly caused by the variation of r0.9. To quantify the spatial extent of the mixing region, the mixed thickness, defined as the distances between the CO 2 fraction level of 10% and 90% along the diagonal line crossing both the CO 2 injection well and the CH 4 production well (see in Figure 7), were compared in Figure 8 under various initial concentrations of CO 2 the cushion gas. The mixed thickness was used to characterize the interface between the working and cushion gases [19]. Clearly, the thickness of the mixed zone increased with the increasing of the CO 2 initiated concentration. The curve of the mixed zone's thickness for the UGSR with initial CO 2 concentration ranging from 8% to 10% was similar, in which the representative curve for the initial 9% of CO 2 could be divided into four stages over the lifetime of project (see in Figure 9c). Besides, the distance between the CO 2 injection well and 10% CO 2 (r 0.1 ), a 90% CO 2 (r 0.9 ) concentration point along the same diagonal line is plotted in Figure 9a,b against time, respectively. The variation of r 0.1 is harmonized with the periodical operation scenarios (Figure 3). It shows that r 0.1 was mainly driven by the gas displacement effect. However, this value was limited in a narrow range Energies 2020, 13, 576 9 of 20 from 189.9 to 193.7 m. Therefore, it could be concluded that the abrupt change of mixed zone thickness (r 0.1 − r 0.9 ) was mostly caused by the variation of r 0.9 . The stage I corresponds to the first injection in the first cycle. In which, r0.1 increased gradually (Figure 9a), while r0.9 fell very fast from 123.5 to 64.4 m (Figure 9b), thus the thickness of the mixed zone increased dramatically due to the mixing of CO2 and CH4 driven by the pressure gradient and the diffusion effect. In stage II, r0.9 decreased only by 4.5 m and r0.1 slightly decreased by about 2 m, so the thickness changes in this stage was only 2.5 m in 245 days, which corresponds to the period from the first well shutting to the second well shutting stage in the first operation cycle. Therefore, the thickness in the second stage behaved smoothly. Although the change of r0.9 was not significant, the concentration of CO2 near the CO2 injection well decreased gradually. With the continuous mixing of CO2 and CH4, the concentration of CO2 in the UGSR decreased to a value lower than 90%, thereby r0.9 decreased to 0 and the distance of the mixed zone substantially increased again in a short time, which was assigned to stage III. In stage IV, the distance of the mixed zone equals r0.1 and changes cyclically with the injection and production of CH4.  The stage I corresponds to the first injection in the first cycle. In which, r0.1 increased gradually (Figure 9a), while r0.9 fell very fast from 123.5 to 64.4 m (Figure 9b), thus the thickness of the mixed zone increased dramatically due to the mixing of CO2 and CH4 driven by the pressure gradient and the diffusion effect. In stage II, r0.9 decreased only by 4.5 m and r0.1 slightly decreased by about 2 m, so the thickness changes in this stage was only 2.5 m in 245 days, which corresponds to the period from the first well shutting to the second well shutting stage in the first operation cycle. Therefore, the thickness in the second stage behaved smoothly. Although the change of r0.9 was not significant, the concentration of CO2 near the CO2 injection well decreased gradually. With the continuous mixing of CO2 and CH4, the concentration of CO2 in the UGSR decreased to a value lower than 90%, thereby r0.9 decreased to 0 and the distance of the mixed zone substantially increased again in a short time, which was assigned to stage III. In stage IV, the distance of the mixed zone equals r0.1 and changes cyclically with the injection and production of CH4. The stage I corresponds to the first injection in the first cycle. In which, r 0.1 increased gradually (Figure 9a), while r 0.9 fell very fast from 123.5 to 64.4 m (Figure 9b), thus the thickness of the mixed zone increased dramatically due to the mixing of CO 2 and CH 4 driven by the pressure gradient and the diffusion effect. In stage II, r 0.9 decreased only by 4.5 m and r 0.1 slightly decreased by about 2 m, so the thickness changes in this stage was only 2.5 m in 245 days, which corresponds to the period from the first well shutting to the second well shutting stage in the first operation cycle. Therefore, the thickness in the second stage behaved smoothly. Although the change of r 0.9 was not significant, the concentration of CO 2 near the CO 2 injection well decreased gradually. With the continuous mixing of CO 2 and CH 4 , the concentration of CO 2 in the UGSR decreased to a value lower than 90%, thereby r 0.9 decreased to 0 and the distance of the mixed zone substantially increased again in a short time, which was assigned to stage III. In stage IV, the distance of the mixed zone equals r 0.1 and changes cyclically with the injection and production of CH 4 .
Unlike the above curve (Figure 9b), there are two stable sections in the r 0.9 curve for the UGSR with 20% CO 2 (Figure 10), which corresponds to the two stable sections of the curve of the mixed zone ( Figure 8). Figure 11 shows the CO 2 concentration in the diagonal of the reservoir for the UGSR with 20% CO 2 . It can be seen that the mixed zone thickness maintained at 131.1 m from 0.35 to 0.84 a, corresponding to the first stable section in Figure 8. During this period, the concentration of CO 2 in the region away from the CO 2 injection well decreased to the threshold value of 90% gradually, then led to an increasing of the mixed zone. It also shows that the mixed zone thickness changed only from 189.7 to 195.9 m during the period of 0.95 to 2.22 a, corresponding to the second stable section in Figure 8.
During this period, the concentration of CO 2 in the region close to CO 2 injection well decreased to lower than 90% gradually, then led to the increasing of the mixed zone again. the region away from the CO2 injection well decreased to the threshold value of 90% gradually, then led to an increasing of the mixed zone. It also shows that the mixed zone thickness changed only from 189.7 to 195.9 m during the period of 0.95 to 2.22 a, corresponding to the second stable section in Figure 8. During this period, the concentration of CO2 in the region close to CO2 injection well decreased to lower than 90% gradually, then led to the increasing of the mixed zone again.  In every production cycle, the CO2 concentration in the produced gas mixture (CO2-CH4) over 15 years is shown Figure 12. Obviously, the CO2 concentration in the produced gas mixture increases in every production cycle and reached up to the maximum at the end of each production. The maximum CO2 concentration increased sharply, when CO2 migrated to the production well in the first few years. Then the maximum concentration of CO2 decreased gradually. This is a result of the decreasing amount of CO2 in the UGSR. The maximum amount of the produced CO2 concentration led to an increasing of the mixed zone. It also shows that the mixed zone thickness changed only from 189.7 to 195.9 m during the period of 0.95 to 2.22 a, corresponding to the second stable section in Figure 8. During this period, the concentration of CO2 in the region close to CO2 injection well decreased to lower than 90% gradually, then led to the increasing of the mixed zone again.  In every production cycle, the CO2 concentration in the produced gas mixture (CO2-CH4) over 15 years is shown Figure 12. Obviously, the CO2 concentration in the produced gas mixture increases in every production cycle and reached up to the maximum at the end of each production. The maximum CO2 concentration increased sharply, when CO2 migrated to the production well in the first few years. Then the maximum concentration of CO2 decreased gradually. This is a result of the decreasing amount of CO2 in the UGSR. The maximum amount of the produced CO2 concentration In every production cycle, the CO 2 concentration in the produced gas mixture (CO 2 -CH 4 ) over 15 years is shown Figure 12. Obviously, the CO 2 concentration in the produced gas mixture increases in every production cycle and reached up to the maximum at the end of each production. The maximum CO 2 concentration increased sharply, when CO 2 migrated to the production well in the first few years. Then the maximum concentration of CO 2 decreased gradually. This is a result of the decreasing amount of CO 2 in the UGSR. The maximum amount of the produced CO 2 concentration in the whole operation process was found at an earlier time with the increasing of the CO 2 cushion gas. It occurred at the 6th and 10th year for the UGSR with 20% and 9% CO 2 cushion gas, respectively. According to the national natural gas standards of China, the concentration of CO 2 in the second class natural gas for civil fuel is not allowed higher than 3% [28]. Therefore, it can be found in Figure 12, that the optimal CO 2 concentration for the cushion gas in this operation scenario was 9%, which would be used as the base case in following studies.
gas. It occurred at the 6th and 10th year for the UGSR with 20% and 9% CO2 cushion gas, respectively. According to the national natural gas standards of China, the concentration of CO2 in the second class natural gas for civil fuel is not allowed higher than 3% [28]. Therefore, it can be found in Figure 12, that the optimal CO2 concentration for the cushion gas in this operation scenario was 9%, which would be used as the base case in following studies.  Temperature   Figures 13-15 show the average reservoir pressure, mixed zone thickness, and CO2 concentration in produced gas over 15 years of operation under different temperatures, respectively. Logically, temperature had a positive influence on reservoir pressure change, and the cyclical changes were synchronized with the operation scenarios ( Figure 3). Figure 14 shows that the thickness of the mixed zone decreased with the rising temperature in the initial stage, which is in accordance with the results from Ma et al. [17]. This result states that the high reservoir pressure caused by increasing temperature increased the dynamic viscosity of the gases. Likewise, the rapidly increased stage (stage III in Figure 9) occurred earlier with a lower temperature. After that, the ranking of the mixed zone thickness reversed eventually in stage IV. The mixing region increased with increasing temperature due to the high diffusion coefficient. Figure 15 shows that the impacts of temperature on the produced CO2 concentration. It can be seen that the effect of temperature on the CO2 concentration in the produced gas could be neglectable ranging from 30 to 50 °C.  Temperature   Figures 13-15 show the average reservoir pressure, mixed zone thickness, and CO 2 concentration in produced gas over 15 years of operation under different temperatures, respectively. Logically, temperature had a positive influence on reservoir pressure change, and the cyclical changes were synchronized with the operation scenarios ( Figure 3). Figure 14 shows that the thickness of the mixed zone decreased with the rising temperature in the initial stage, which is in accordance with the results from Ma et al. [17]. This result states that the high reservoir pressure caused by increasing temperature increased the dynamic viscosity of the gases. Likewise, the rapidly increased stage (stage III in Figure 9) occurred earlier with a lower temperature. After that, the ranking of the mixed zone thickness reversed eventually in stage IV. The mixing region increased with increasing temperature due to the high diffusion coefficient. Figure 15 shows that the impacts of temperature on the produced CO 2 concentration. It can be seen that the effect of temperature on the CO 2 concentration in the produced gas could be neglectable ranging from 30 to 50 • C.

Effect of the Reservoir
Energies 2020, 13, 576 13 of 21 . Figure 13. Reservoir average pressure for different temperature.     . Figure 13. Reservoir average pressure for different temperature.

Effect of the Reservoir Permeability
The permeability of depleted gas reservoirs may range from a few to hundreds of millidarcy. To quantify the impact of reservoir permeability on the mixing behavior of the gases, some cases with the reservoir horizontal permeability of 50, 70, 100, and 120 mD were carried out. As shown in Figure 16, the reservoir permeability had only a minor effect on the formation pressure, but the effect on the mixing behavior of CO 2 and CH 4 was much more significant. Figure 17 shows that the mixing thickness was positively correlated to permeability and this relationship was inversed after entering stage IV. In the production stage of the first year (stage I in Figure 9c), higher permeability will accelerate the migration of CO 2 , further promoting the mixing behavior. With continuous operation, more CO 2 was produced along with CH 4 production in the case with higher permeability, which made the CO 2 concentration point of 10% quickly return to the CO 2 injection well. In other words, r 0.1 and the corresponding thickness of the mixed zone decreased much more quickly with higher permeability. Similar to that mentioned above, the thickness of the mixed zone changed periodically. The gases were much easier to mix in the reservoir with higher permeability. In stage IV (Figure 9c), the mixing zone was dominantly determined by r 0.1 , the displacing effect of CH 4 injection was more significant in the permeable reservoir, thus the mixing zone was decreased with increasing permeability and varied in a large amplitude. permeability of 100 and 120 mD affected the distribution of CO2 in the reservoir a lot and led to the non-uniform fluctuation of mixed region as shown in Figure 17. Therefore, the CO2 concentration in a highly permeable reservoir might be less than that of a lowly permeable reservoir in the late period of an operation if no CO2 was re-injected for supplementation. Likewise, the maximum CO2 concentration was found at the earlier time in a highly permeable reservoir.    Figure 19 depicts the formation pressure for the reservoirs with different residual water saturations. The formation pressure increased with the growing residual water saturation, because The CO 2 concentration in the produced gas for different reservoir horizontal permeability was compared in Figure 18. Clearly, CO 2 concentration in produced gas was positively correlated to the permeability. The high CO 2 production rates between 3 and 8 years for the reservoir with a permeability of 100 and 120 mD affected the distribution of CO 2 in the reservoir a lot and led to the non-uniform fluctuation of mixed region as shown in Figure 17. Therefore, the CO 2 concentration in a highly permeable reservoir might be less than that of a lowly permeable reservoir in the late period of an operation if no CO 2 was re-injected for supplementation. Likewise, the maximum CO 2 concentration was found at the earlier time in a highly permeable reservoir.   Figure 19 depicts the formation pressure for the reservoirs with different residual water saturations. The formation pressure increased with the growing residual water saturation, because uncompressible water can build up higher pressure in the reservoir. Figure 20 shows the thickness of the mixed zone. In the production stage of the first year (stage I in Figure 9c), the mixing zone in the reservoir tended to increase with the expanding residual water saturation, one reason is that the occupied pores promote the injected CO2 traveling deeper into the reservoir, i.e., higher mixing thickness under the same injection volume. Another reason is that the residual water also generates more tortuous flow-paths for the gases as shown in Figure 21, which was discovered in the result of nuclear magnetic resonance during a core flood experiment [29]. Therefore, the rapidly increasing stage (stage III in Figure 9c) occurred earlier with higher residual water saturation. For the same reason, the mixing thickness increasing happened firstly, in the case with higher residual water saturation after initiation.  Figure 19 depicts the formation pressure for the reservoirs with different residual water saturations. The formation pressure increased with the growing residual water saturation, because uncompressible water can build up higher pressure in the reservoir. Figure 20 shows the thickness of the mixed zone. In the production stage of the first year (stage I in Figure 9c), the mixing zone in the reservoir tended to increase with the expanding residual water saturation, one reason is that the occupied pores promote the injected CO 2 traveling deeper into the reservoir, i.e., higher mixing thickness under the same injection volume. Another reason is that the residual water also generates more tortuous flow-paths for the gases as shown in Figure 21, which was discovered in the result of nuclear magnetic resonance during a core flood experiment [29]. Therefore, the rapidly increasing stage (stage III in Figure 9c) occurred earlier with higher residual water saturation. For the same reason, the mixing thickness increasing happened firstly, in the case with higher residual water saturation after initiation.  Like the effect of permeability, the mixing thickness was inversely correlated to the residual water saturation after entering stage IV. This was caused by the dissolution of CO 2 in residual water. As shown in Figure 22, less CO 2 existed as gas in the reservoir with higher water saturation, due to more CO 2 dissolved in the residual water. Figure 22 also shows the dissolved CO 2 changed periodicity in a certain range during operation. These are results from the dissolution of CO 2 , which is principally pressure-and temperature-dependent. The CO 2 concentration in the produced gas for different residual water saturation is shown in Figure 23. It demonstrated that the effect of the residual water on the CO 2 concentration in the produced gas was still minor, which was due to the differences of the gaseous CO 2 in the UGSR being not very significant as shown in Figure 22.   Like the effect of permeability, the mixing thickness was inversely correlated to the residual water saturation after entering stage IV. This was caused by the dissolution of CO2 in residual water. As shown in Figure 22, less CO2 existed as gas in the reservoir with higher water saturation, due to more CO2 dissolved in the residual water. Figure 22 also shows the dissolved CO2 changed    Like the effect of permeability, the mixing thickness was inversely correlated to the residual water saturation after entering stage IV. This was caused by the dissolution of CO2 in residual water. As shown in Figure 22, less CO2 existed as gas in the reservoir with higher water saturation, due to more CO2 dissolved in the residual water. Figure 22 also shows the dissolved CO2 changed  [29].

Effect of the Residual Water
Energies 2020, 13, 576 17 of 21 periodicity in a certain range during operation. These are results from the dissolution of CO2, which is principally pressure-and temperature-dependent. The CO2 concentration in the produced gas for different residual water saturation is shown in Figure 23. It demonstrated that the effect of the residual water on the CO2 concentration in the produced gas was still minor, which was due to the differences of the gaseous CO2 in the UGSR being not very significant as shown in Figure 22.

Effect of the Production Rate
The average reservoir pressure under different production rates were plotted against time in Figure 24. Noticeably, average reservoir pressure changed periodically and decreased with the increasing production rate. Comparing the mixing thickness shown in Figure 25, the impact on the mixing behavior of CO2 and CH4 was limited, especially in the early operation stages I-III. Considering the relatively significant difference of the mixed zone at the same operation rate (same pressure gradient) while at a different temperature and residual water saturation (different diffusion effect) in Figures 14 and 20, it could be inferred that the mixing of the gases in these stages is controlled by the diffusion effect rather than pressure gradient. In the subsequent operation, stage IV, the mixing zone had a slight decrease with an increasing production rate due to the impact of the

Effect of the Production Rate
The average reservoir pressure under different production rates were plotted against time in Figure 24. Noticeably, average reservoir pressure changed periodically and decreased with the increasing production rate. Comparing the mixing thickness shown in Figure 25, the impact on the mixing behavior of CO 2 and CH 4 was limited, especially in the early operation stages I-III. Considering the relatively significant difference of the mixed zone at the same operation rate (same pressure gradient) while at a different temperature and residual water saturation (different diffusion effect) in Figures 14 and 20, it could be inferred that the mixing of the gases in these stages is controlled by the diffusion effect rather than pressure gradient. In the subsequent operation, stage IV, the mixing zone had a slight decrease with an increasing production rate due to the impact of the pressure gradient. As shown in Figure 26, production also has a positive influence on the produced CO 2 . The maximum amount of produced CO 2 was limited to under 3%, which still satisfies the required standard [28]. Therefore, the production rate at 4.2 kg/s could be used in this case, to improve the efficiency of the UGSR.
Energies 2020, 13, 576 18 of 21 pressure gradient. As shown in Figure 26, production also has a positive influence on the produced CO2. The maximum amount of produced CO2 was limited to under 3%, which still satisfies the required standard [28]. Therefore, the production rate at 4.2 kg/s could be used in this case, to improve the efficiency of the UGSR.

Conclusions
The Donghae gas reservoir was used to study the suitability of utilizing CO2 as the cushion gas for a depleted gas reservoir transformed gas storage reservoir. The maximum amount of CO2 concentration for cushion gas was 9% under the condition of production and injection for 120 d and 180 d in a production cycle at a rate of 4.05 kg/s and 2.7 kg/s, respectively.
The typical curve of the mixing zone thickness could be divided into four stages. They were the increasing stage, the smooth stage, the suddenly increasing stage, and the periodic change stage. It should be mentioned that there existed two smooth stages and two suddenly increasing stages when the CO2 concentration in the cushion gas increased to 20%. In the periodic change stage, the mixed zone thickness increased with the increasing of the CO2 fraction, temperature, production rate, and the decreasing of permeability and water saturation. The correlation of the mixing zone thickness with reservoir temperature, permeability, and residual water was inverse in the former stages.
The CO2 fraction in the cushion gas, reservoir permeability, and production rate had a significant effect on the breakthrough of CO2 in the production well, while the effect of water saturation and temperature was relatively minor. For the purpose of utilizing more CO2 as cushion gas in the UGSR, CO2 should be injected for supplementation during the operation of UGSR, which will be discussed in our future work.
Author Contributions: C.C. and Z.H. proposed the main idea; C.C. and J.L. conducted the simulation and performed the analysis; C.C. and J.L. wrote the paper; H.X., F.M. and X.W. contributed to the development and discussion of the results.

Conclusions
The Donghae gas reservoir was used to study the suitability of utilizing CO 2 as the cushion gas for a depleted gas reservoir transformed gas storage reservoir. The maximum amount of CO 2 concentration for cushion gas was 9% under the condition of production and injection for 120 d and 180 d in a production cycle at a rate of 4.05 kg/s and 2.7 kg/s, respectively.
The typical curve of the mixing zone thickness could be divided into four stages. They were the increasing stage, the smooth stage, the suddenly increasing stage, and the periodic change stage. It should be mentioned that there existed two smooth stages and two suddenly increasing stages when the CO 2 concentration in the cushion gas increased to 20%. In the periodic change stage, the mixed zone thickness increased with the increasing of the CO 2 fraction, temperature, production rate, and the decreasing of permeability and water saturation. The correlation of the mixing zone thickness with reservoir temperature, permeability, and residual water was inverse in the former stages.
The CO 2 fraction in the cushion gas, reservoir permeability, and production rate had a significant effect on the breakthrough of CO 2 in the production well, while the effect of water saturation and temperature was relatively minor. For the purpose of utilizing more CO 2 as cushion gas in the UGSR, CO 2 should be injected for supplementation during the operation of UGSR, which will be discussed in our future work.