Chuan-Kai Niu，Yu-Fei Tan

(School of Municipal＆Environment Engineering，Harbin Institute of Technology，Harbin 150090，China)

1 Introduction

The technology of using CO2as cushion gas for UGSR has many advantages，including replacing the large deposits of natural gas in UGSR，achieving carbon sequestration，reducing the greenhouse effect，and bringing high economic benefits，etc［1］.According to the statistics from United States in 2010［2-4］，the total capacity of cushion gas in its running 411 UGSR amounts to 117.67×109m3.If 20% of the cushion gases are replaced with CO2，as much as about 23.54× 109m3natural gas deposits can be further utilized，which bringsthe significantsocialand economic benefits［5］.

The major technical challenge of using CO2as the cushion gas in reservoirs is the mixed of CO2and natural gas.The studies have shown that the mixed of two gases will cause the increase of impurities in the mined peak shaving gas and the decrease of the calorific value［6-7］.The mixed zone of CO2and natural gas wascaused by both moleculardiffusion and convection flow.However，the previous mathematic models［8-9］often neglected the influence of molecular diffusion.In these models，the coefficient for the molecular diffusion term is much smaller than the coefficient for the convection diffusion term.So，these models cannot accurately predict the mixing extent of the two gases and the change of these concentrations.

In this paper，a numerical model is proposed to predict the mixture and diffusion of CO2and natural gas in aquifer UGSR，based on the three dimensional gaswater two-phase flow theory and the gas diffusion theory. Firstly， the pressure equation and the saturation equation for the reservoir are discretized and solved using finite difference method，so the transient pressure distribution and saturation in the reservoir are obtained.Then implicit difference scheme is used to discretize the gas diffusion equation，through which the concentration distribution of CO2and natural gas of reservoir are derived based on the transient pressure distribution of reservoir and saturation.The accuracy of the model is further validated using the actual operation data from Dazhangtuo UGSR，China.Using the validated model，the mixed characteristic of CO2and natural gas and the migration mechanism of the mixed zone in an underground porous reservoir is further studied. Particularly，the impacts of the following factors on the migration mechanism are studied:the ratio of CO2injection，the reservoir porosity and the initial operation pressure.Based on the results，the optimal CO2injection ratio and an optimal control strategy to manage the migration of the mixed zone are obtained. These results will provide the technical guides for using CO2as cushion gas for UGSR in real projects.

2 Development and Solution of Mathematical Model

2.1 Gas-Water Two-Phase Seepage Equation

Following assumptions have been made to simplify the mathematical model:(1)temperature is uniformly distributed in the horizontal direction;(2)natural gas is insoluble.Since the seepage is governed by gravity and capillary pressure and italso considersthe compressibility of rock and fluid，and the anisotropy of the rock，theseepagecan bedescribed bythe extended Darcy Law［10-11］.

Gas phase:

Water phase:

where K is absolute permeability of mixed gas(μm2); Krgand Krware relative permeability of mixed gas and water phase;μgand μware dynamic viscosity of mixed gas and water phase(Pa·s);Sgand Sware gas and water saturation;φ is effective porosity;δ is a function of gas source and δ=1 at the wellhead and δ=0 not at the wellhead;H is dimension factor;Cgand Cware isothermal compressibility of mixed gas and water phase，which can be calculated as in Eq.(3)，and Bgand Bware volume coefficients of mixed gas and water phase，which can be calculated as in Eq.(4).

where Vgand ρgare the volume(m3)and density of mixed gas(kg/m3)under formation condition;Vgsand ρgsarethevolume(m3)and densityofmixed gas(kg/m3)under standard condition;Vwand ρware the volume(m3)and density of water(kg/m3)under formation condition;Vwsand ρwsare the volume(m3) and density ofwater(kg/m3)understandard condition;pgand pware capillary pressure of mixed gas and water phase(Pa);qgand qware unit injectionproduction rate of mixed gas and water phase at wellhead(kg/(m3·s))，and it is positive while injection and negative while production and zero while the well is shut-in.

2.2 Diffusion Equations of Mixed Gas

The underground reservoir has good fractal characters and the mixed gas diffusion follows Fick＇s law.The total mass flow rate of the gas diffusion process contains convection and molecular diffusion. Convective flow of fluid is the mass transfer caused by seepage flow of mixed gas，while molecular diffusion is the molecular motion from high concentration range to low concentration range caused by concentration difference of mixed gas.Considering the effects of dimension factor and coupling mass， diffusion equations of mixed gas are expressed in Eqs.(5)and (6)［12-13］，in which A and B represent CO2and natural gas respectively.

where DAB，ris fractal diffusion coefficient of CO2in the medium of natural gas(m2/s)，while，DBA，ris fractal diffusion coefficient of natural gas the medium of CO2(m2/s);MAand MBare molecular weight of CO2and natural gas;CAand CBare molar concentration of CO2and natural gas(mol/m3).

2.3 Auxiliary Equations

Saturation equilibrium equation:

The constraint equation of capillary pressure which is a function of fluid saturation is as follows:

Now，it can complete the variable resolution of pg，pw，Sg，Sw，CAand CBusing the coupled equations of Eqs.(1，2，5-8).Moreover，the coefficient of each parameter in the equations needs to be calculated based on the auxiliary equations，in which the detailed derivation processes are shown in Ref.［14］.

2.4 Initial Conditions and Boundary Conditions

The first boundary condition is used at the outer boundary of UGSR as follows:

However，the second boundary condition is given at the head of the injection-production well of UGSR as follows，in which n is normal direction.

The boundary condition of the UGSR in the sandwich is given as follows:

The initial pressure and initial saturation of UGSR which are certain functions are given as follows:

The initial pressure gradient in the reservoir of the UGSR is given as follows:

2.5 Solution of Numerical Model

During the operation of UGSRs in aquifer with cushion gas，the edge and bottom water are not usually to contact with the mixed zone.Therefore，the pressure fields and concentration fields of natural gas and cushion gas are respectively solved in the pure natural gas field and pure cushion gas field，while the coupled pressure fields and concentration fields in the mixture gas field are solved using the iterative method.

Inserting saturation equilibrium Eq.(7)and capillary pressure Eq.(8)into gas-water two-phase seepage Eqs.(1-2)and then eliminating the variables pwand Sw，the simultaneous equations of pgand Sgare obtained as follows:

The reservoir space of UGSR is a type of porous media with fractal characteristics［15］. The relative permeability of gas phase which are not fixed values in the above equations are changing with the difference of pore fractal dimension and porosity of porous media. The relative permeability of gas phase Krgis calculated as in Eq.(17)［16-17］.

where dfis pore fractal dimension of porous media;rminand rmaxare minimum pore size and maximum pore size (μm);λ is aperture ratio of porous media，and λ= rmin/rmax;C is a constant related with the structure of reservoir space.

Reservoir pressure equation and saturation equation are discretized using finite element method，and the discrete form of the equations is as follows:

where［K］，［R］and［E］are rigidity matrix;［G］，［N］，［F］and［O］are coefficient matrix.

Based on the calculation method of the binary spread system，gas diffusion equations are discretized and solved using finite difference method，where，the fractal diffusion coefficient in porous media(Di，r)is calculated as in Eq.(19)［18］.

where Di，ais diffusion coefficient of CO2or natural gas in the homogeneous Euclidean space(m2/s)，and i represent CO2and natural gas respectively;θ is diffusion attenuation index and θ=1/H-2，here H called Hurst coefficient is a characteristic parameter reflected diffusion slowing-down effect of molecular Brownian motion in different fractal structure［19］.

Taking CO2for example，first of all，Eq.(5)can be converted into the following form:

Non-isometric 8 points implicit scheme of Eq. (20)is discretized using finite difference method as follows.Where，τ is time step.

Similarly，the diffusion equation of natural gas is discretized. Firstly， it is assumed the initial concentration of reservoir，the pressure field and saturation field at the next time step are calculated according to the equations of the discrete pressure and saturation by finite element method. Then the concentration field at this time step is calculated using the discrete diffusion equation based on the calculation result of pressure field and saturation field. The calculation of pressure field and saturation field at this moment is based on the past moment concentration field，and vice versa.In addition，the calculation of characteristic parameters is based on the result of pressure field，saturation field and concentration field. The detailed calculation steps to solve the coupling of pressure field，saturation field and concentration field are as follows:

1)To give the calculated value or initial value (Pg，i，yi)of node pressure and mole fractions of working gas and cushion gas(hereinafter short for twoparameter)and the related formation parameter.

2)According the calculated value or initial value (Pg，i，yi)，characteristic parameter obtains by chart method.The pressure field and the saturation field of next momentand(iteration expression with a superscript)are calculated by pressure equation (15)and(16).The velocity distributionof the pressure according to the calculated value or initial value Pg，i.And the concentration distribution，)are calculated byEq.(21)and initial concentration(CA，i+1，CB，i+1)by the calculated value or initial value(Pg，i，yi).Then the mole fractions of next timeof the two-parameter.

3)Taking the calculated value or initial valueand the calculated valueof Step 2 on average，denotedrepeating Step 2 is taking average as the condition of the next iterative calculations and obtained the calculated value

4)Repeating Step 3，but taking the calculated valuewhich is the average ofandand initial concentration distribution (CA，i+1，CB，i+1)as the input conditions.The iterations are completed once meeting the conditions of|and，and the calculated valueis the input condition of next step.

5)Taking the result of two-parameterdenote)，which is the input condition of the time step according the steps from 1 to 4.

3 Model Validation

Fig.1 shows the well-site distribution and reservoir sketch of Dazhangtuo UGSR in Tianjin，China.The well marked as‘tuozhu-2’is the central injection-production wellaccording to the actual operation and the wellarray maps［20］.The original formation pressure and temperature are 24.17 MPa and 378.15 K.According to the high pressure mercury injection experimentfor a rock sample ofthis reservoir［21-22］，the convective diffusion coefficient of the rock in the homogeneous Euclid space is Da= 0.517 cm2/s.The average permeability of rock sample of Dazhangtuo reservoir at several different wells is K= 15.31×10-3μm2，according to the core permeability measured data of different reservoir of Dazhangtuo underground gas storage［23］. The pore structure parameters of the rock sample and other physical parameters used in the calculations are shown in Table 1.

Fig.1 Well-site distribution and reservoir sketch of Dazhangtuo UGSR

Table 1 Fractal structure parameters of rock swatch and other properties of the model in UGSR

Fig.2 shows the real measurement data［20，23］and the simulation results of daily gas production volume and reservoir pressure of the‘tuozhu-2’well in two periods from November 15th，2000 to March 15th， 2001 and from November 15th，2001 to March 15th，2002.And there is no cushion gas injection in the first period and there is in the second period.In the simulation，the time step is 0.5 days，and the space is divided into small grids with width，length and height of 1 m，1 m and 0.1 m respectively.

Fig.2 Changes of gas production and reservoir pressure with or without CO2injection from edge well

The comparison of measured and simulated daily produced gas volume in the‘tuozhu-2’well in two gas production periods is shown in Fig.2(a).It shows that in the case without cushion gas，the maximum absolute error of daily produced gas volume between measurement and simulation is 4.2×103m3.In the case with cushion gas injection，the maximum absolute error of daily produced gas volume is 8.3×103m3.As shown in Fig.2(a)，the daily produced gas volume is little difference between simulation predicted results and actual operation data，which can meet the project requirement.Fig.2(b)shows the reservoir pressure in two gas production periods.The maximum relative error of reservoir pressure drop is 13.2% in the gas production cycle without cushion gas injection.And the maximum relative error of reservoir pressure drop is 14.3%in the gas production cycle with cushion gas injection.The pressure drop in the simulation is always lower than that in actual production operation process because of the change of CO2solubility in the aquifer. Therefore，CO2solubility under instantaneous pressure should be iterative corrected by reservoir pressure and temperature in the simulation，in order to reduce the error between the simulation predicted results and actual operation data.The above comparison shows that the gasmixture prediction modelhasgood accuracy.Therefore，the developed model can be used to simulate the operation process of the UGSR with cushion gas.

4 Simulation Results and Sensitivity Analysis of Mixed Zone

To perform the sensitivity analysis，a simplified version of the Dazhangtuo UGSR is used.In this model，the‘tuozhu-2’well is located in the center， the‘tuozhu-1’well，the‘ban53’well，the‘ban57’well and the‘banshen3-1’well are distributed evenly around it.The area of the UGSR in this model is 2 km wide by 2 km，and the height is 20 m.Daily injected or produced gas volume of the central well is 350×103m3.The original formation pressure and temperature are 24.17 MPa and 378.15 K，and all of the boundaries are assumed to be closed. In the simulation，the time step is 0.5 days，and the space is divided into small grids with width，length and height of 1 m，1 m and 0.1 m respectively.The schedule for a complete injection-production cycle is shown in Table 2.Other parameters used in the analysis are shown in Table 2.

Table 2 Working process of a complete injectionproduction cycle

4.1 Ratio of CO2Injection in Cushion Gas

Fig.3 shows the impact of the CO2ratio to the change ofthe mixed zone atdifferentinjectionproduction cycles.At the injection initial phase，when a certain amount of CO2has been detected near the center well，the formation of the mixed zone has been accelerated.At the end of the fourth gas production cycle，the volume ratio of impurities(CO2)in the mined the peak shaving gas from production well has reached 10.4%，which results in the decrease of gas calorific value and cannot satisfy the requirement of the real project.When 30% CO2(volume ratio)is injected at the edge well at the initial injection phase，the mixing process of natural gas and CO2is slowed. At the end of the fourth gas production cycle，CO2has been not detected in the production gas in the central well.This shows that injecting CO2cushion gas from the edge wells can effectively slow down the process of the formation and migration of mixed zone，and thus is suggested as the injection method to be used in practice.

Fig.3 Mixed zone in different injection-production cycles

Fig.4 shows that when injecting 30% CO2(volume ratio)as cushion gas at the initial injection phase，CO2has been detected at the bottom of the injection-production well after the tenth production cycle.However，when all the gas in reservoir is CO2(100%)at the initial injection phase，CO2can be detected at the bottom of the injection-production well as early as the end of the third production cycle.This seriously affects the quality of the peak shaving gas from production well.Fig.5 shows the impact of the different initial CO2injection ratio at the edge wells to the width of the natural gas-CO2mixed zone after one production cycle.As the CO2injection ratio decreases，the width of the mixed zone also decreases.After the CO2injection ratio decreases from 90%to 30%，the width of mixed zone decreases from 271.8 m to 134.4 m.Balancing the operation stability of the UGSR and the operation cost，the optimal CO2injection ratio at the edge well is 25.6%.

Fig.4 Cycle that cushion gas always appears at the bottom of I/P well

Fig.5 Width of mixed zone under different CO2initial ratio of cushion gas

4.2 Reservoir Porosity

Fig.6 shows the width and the location of the mixed zone at different porosities when 30% CO2(volume ratio)is injected from the edge well at the initial injection phase.As shown in Fig.6，as the porosity increases，the width ofthe mixed zone decreases，and the distance to the centralwell increases.This result suggests that large porosity can slow down the mixing of natural gas and CO2and provide larger space for storage.Therefore，reservoir with large porosity is better suited as UGSR.Fig.6(b) shows that the mixing extent of natural gas and CO2is high at the initial injection phase.As natural gas is injected，the width of the mixed zone decreases and after the injection phase is finished，the width of the mixed zone gradually increases and the growth has slowed.At the initial injection phase，the mixing between natural gas and CO2is quick and the width grows rapidly.However，as the volume of natural gas increases，a pure natural gas region is formed around the central well，which drives the mixed zone to move outward away from the central well.Thus，the width of the mixed zone decreases.When the injection phase is finished，the mixed zone develops stably and the width gradually increases again.

Fig.6 Migration of mixed zone under different porosities

4.3 Initial Operation Pressure

Fig.7 shows the impact of the initial pressure on the width and location of the mixed zone.It is seen from Fig.7 that，as the initial pressure decreases，the mixing extent of natural gas and CO2increases，After the initial pressure decreases from 16 MPa to 12 MPa，the distance of the mixed zone to the central well after an operation period of 220 days increases from 217.5 m to 281.4 m(increases by 29.4%)，and the average width increases from 129.5 m to 154.9 m (increases by 19.6%).Thus，the initial pressure of the reservoir should be as large as possible.However，it should be noted that as the initial pressure is related with the ratio of the working volume to the total volume，the initial pressure is bounded by a reasonable region.Usually，the initial pressure is slightly larger than the minimally allowed pressure in the reservoir.

Fig.7 Migration of mixed zone under different Initial pressure

5 Conclusions

A numerical model based on the theory of three dimensional gas-water two-phase flows and the theory of gas diffusion is developed in this paper to study the usage of CO2as cushion gas in UGSR.It is found that the model produces simulation results close to the operation data of a real UGSR.Using this model，the growth and change of the natural gas-CO2mixed zone is studied，and the following observations are drawn:

1)When all of the cushion gas in the reservoir is CO2，the mixed zone of natural gas and CO2will form rapidly near the injection-production well and the operation of the UGSR will be seriously affected.When CO2is used as only part of the cushion gas and injected into the UGSR from the edge wells，the impact of the mixed zone to the operation will be negligible.A calculation through multi-cycle injection-production simulation shows that the optimal CO2injection ratio is 25.6%.

2)In a reservoir of high porosity，the formation and migration of the mixed zone is slower during its operation.Thus，a reservoir of high porosity is better suited for storing natural gas.

3)A large initial injection pressure helps limiting the formation and migration ofthe mixed zone. However，as increasing the initial injection pressure decreases the storage space in the reservoir，a proper initial injection pressure is crucial.In general，the initial injection pressure of reservoir is slightly higher than the minimally allowed reservoirpressure in UGSR.

［1］Procesi M，Cantucci B，Buttinelli M，et al.Strategic use of the underground in an energy mix plan:Synergies among CO2，CH4geological storage and geothermalenergy. Latium Region casestudy(CentralItaly). Applied Energy，2013，110(10):104-131.

［2］Wang Baohui，Yan Xiangzhen，Yang Xiujuan.Study on the dynamic migration law of inert cushion gas during injection-production process in underground natural gas storage reservoir.Science Technology and Engineering，2013，13(31):9184-9188.

［3］Yin Huchen，Chen Junbin，Lan Yifei，et al.Technical development status quo and inspiration of typical gas storages in North America. Oil＆ GasStorage and Transportation，2013，32(8):814-817.

［4］ WangXiuli，EconomidesM J. Purposefullybuilt underground natural gas storage.Journal of Natural Gas Science and Engineering，2012，9(11):130-137.

［5］ Plaat H.Underground gas storage:Why and how. Geological Society，London，Special Publications，2009，313:25-37.

［6］Tan Yufei，Cao Lin，Lin Tao.Feasibility analysis about carbon dioxide as cushion gas for natural gas storage.Oil＆Gas Storage and Transportation，2006，25(3):12-14.

［7］Chen Mingjie，Buscheck T A，Wagoner J L，et al. Analysis of fault leakage from Leroy underground natural gas storage facility.Hydrogeology Journal，2013，21(7): 1429-1445.

［8］ Tan Yufei.Technology and Numerical Simulation of Underground Gas Storage Reservoir.Beijing:Petroleum Industry Press，2007.171-178.

［9］Tan Yufei，Lin Tao.Simulation research on carbon dioxide as cushion gas in gas underground reservoirs.Journal of Harbin Institute of Technology(New Series)，2009，16 (1):87-90.

［10］Leontoev N E.Description of weakly compressible fluid flows in porous media for a nonlinear seepage law.Fluid Dynamics，2013，48(3):402-406.

［11］Navas P，Susana Lopez-Querol.Generalized unconfined seepage flow model using displacement based formulation. Engineering Geology，2013，166(8):140-151.

［12］ Lin Tao. Simulation ofInjection-Production and Optimization of Running Control in Gas Reservoirs Using Carbon Dioxide.Harbin:Harbin Institute of Technology，2008.

［13］Vopiěka O，De Angelis M G，Du Naiying，et al.Mixed gas sorption in glassy polymeric membranes:II.CO2/CH4mixtures in a polymer of intrinsic microporosity(PIM-1). Journal of Membrane Science，2014，459(6):264-276.

［14］Li Juanjuan.Simulation Research of Use of Inert Cushion Gas in Underground Gas Storage Reservoir in Aquifer. Harbin:Harbin Institute of Technology，2006.

［15］Yu Benfu，Yan Xiangzhen，Yang Xiujuan.Seepage law analysis of natural gas considering the fractal features of crack reservoirin underground gas storage. Science Technology and Engineering，2014，14(5):54-61.

［16］Guarracino L，R?tting T，Carrera J.A fractal model to describe the evolution of multiphase flow properties during mineral dissolution.Advances in Water Resources，2014，67(5):78-86.

［17］YangYu，Sun Hansen，PengXiaodong，etal. Quantitative study on fractal characteristics of the structure of CBM reservoir.Special Oil＆Gas Reservoirs，2013，20 (1):31-34.

［18］Emera M K，Sarma H K.Genetic algorithm-based correlations offer more reliable prediction of minimum miscibility pressures between reservoir oil and CO2or flue gas.Journal of Canadian Petroleum Technology，2007，46 (8):19-25.

［19］Li Jin，Huang Jianhua.Dynamics of stochastic non-Newtonian fluids driven by fractional Brownian motion with Hurst parameter H∈ (1/4，1/2).Applied Mathematics and Mechanics(English Edition)，2014，34(2):189-208.

［20］Yang Shuhe，He Shumei，Yang Bo，et al.The operation practice and evaluation for Dazhuangtuo underground gas storage.Natural Gas Geoscience，2003，14(5):425-428.

［21］HuangWeihe，Yang Yuling. Reservoirprotection regarding underground gas storage construction based on low-pressure resrrvoir.Natural Gas Industry，2008，28 (4):102-104.

［22］Mao Zhiqiang，Xiao Liang，Wang Zhaonian，et al. Estimation of permeability by integrating nuclear magnetic resonance(NMR)Logs with Mercury Injection Capillary Pressure(MICP)Data in tight gas sands.Applied Magnetic Resonance，2013，44(4):449-468.

［23］Yang Feng，Ning Zhengfu，Kong Detao，et al.Pore structure of shales from high pressure mercury injection and nitrogen adsorption method. NaturalGasGeoscience，2013，24(3):450-455.

［24］Ma Huifang.Numeric Simulation of Dazhangtuo Gas Storage.Beijing:China University of Geosciences，2002.