Hydromechanical Response and Impact of Gas Mixing Behavior in

Jun 24, 2019 - Feasible injection rates and injection schedules can be derived from the ... The amount of cushion gas generally accounts for 40–70% ...
2 downloads 0 Views 7MB Size
Article Cite This: Energy Fuels 2019, 33, 6527−6541

pubs.acs.org/EF

Hydromechanical Response and Impact of Gas Mixing Behavior in Subsurface CH4 Storage with CO2‑Based Cushion Gas Jianli Ma,†,‡,§ Qi Li,*,†,‡ Thomas Kempka,§,∥ and Michael Kühn§,∥ †

Downloaded via GUILFORD COLG on July 25, 2019 at 09:27:29 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.

State Key Laboratory of Geomechanics and Geotechnical Engineering, Institute of Rock and Soil Mechanics, Chinese Academy of Sciences, Wuhan, 430071, China ‡ University of Chinese Academy of Sciences, Beijing 100049, China § Fluid Systems Modelling, GFZ German Research Centre for Geosciences, 14473 Potsdam, Germany ∥ Earth and Environmental Science, University of Potsdam, 14469 Potsdam, Germany ABSTRACT: Power-to-gas (PtG) stores chemical energy by converting excess electrical energy from renewable sources into an energy-dense gas. Due to its higher available capacity compared to surface-based storage technologies, subsurface storage in geological systems is the most promising approach for efficient and economic realization of the PtG system’s storage component. For this purpose, methane (CH4) produced by methanation by means of hydrogen (H2) and carbon dioxide (CO2) is stored in a geological reservoir until required for further use. In this context, CO2 is used as the cushion gas to maintain reservoir pressure and limiting working gas, i.e., (CH4) losses during withdrawal periods. Consequently, mixing of both gases in the reservoir is inevitable. Therefore, it is necessary to minimize the gas mixing region to optimize the efficiency of the PtG system’s storage component. In the present study, the physical properties of CH4, CO2 and their mixtures are reviewed. Then, a multicomponent flow model is implemented and validated against published data. Next, a hydromechanically coupled model is established, considering fluid flow through porous media and effective stresses to investigate the mixing behavior of both gases and the mechanical reservoir stability. The simulation results show that, with increasing reservoir thickness and dip angle, the mixing region is reduced during gas injection if CO2 is employed as the cushion gas. In addition, the degree of mixing is lower at higher temperatures. Feasible injection rates and injection schedules can be derived from the integrated reservoir stability analysis. The methodology developed in the present study allows the determination of optimum strategies for storage reservoir selection and gas injection scheduling by minimizing the gas mixing region.

1. INTRODUCTION Mitigation of global climate change is one of humankind’s major challenges in this century.1,2 The content of the main driver of the anthropogenically driven greenhouse gas effect, carbon dioxide (CO2), has increased by 40% in the atmosphere over the past 250 years.3,4 Several studies show that carbon dioxide capture and storage (CCS) is one of the most promising approaches for greenhouse gas emission mitigation.5−9 Meanwhile, an increasing number of countries have established substantial efforts to develop the renewable energy sector. However, the renewable energy supply fluctuates on different time scales, mainly due to weather conditions. Large amounts of renewable energy are wasted during periods of low power demand.10,11 As a proven energy conversion technology, power-to-gas (PtG) can contribute to effectively store renewable energy, and thus mitigate energy losses.12−15 PtG combined with subsurface gas storage provides an opportunity for CO2 utilization,16 while methane (CH4) produced from excess electricity by methanation can be temporarily stored below ground, until its combustion is required, to provide energy when renewable sources cannot fulfill the demand. Subsurface natural gas storage has been widely applied due to economic reasons and the high storage capacities available in saline aquifers and depleted gas reservoirs.17 PtG, combined with subsurface storage, enables efficiently mastering the divergence between electricity supply © 2019 American Chemical Society

and demand, reducing the waste of excess energy from renewables as well as CO2 emissions.18,19 However, natural gas underground storage requires a cushion gas to maintain the required reservoir pressure, prevent water intrusion, and ensure the stability of gas storage. The amount of cushion gas generally accounts for 40−70% of the total gas storage volume.20 At present, gas storage in most countries worldwide uses natural gas as a cushion gas. When the storage site is abandoned, a considerable amount of cushion gas is trapped, resulting in a large amount of energy loss. For that reason, inert gases have been considered to replace CH4 as a cushion gas, with successful commercial-scale application, as demonstrated in France in the 1980s.20 Since then, several publications have discussed the efficiency of using an inert cushion gas, especially nitrogen (N2).21−24 However, due to the small density difference between N2 and CH4, gas mixing is quite severe during the withdrawal. In addition, a cushion gas with relatively low compressibility will cause a sharp pressure increase during injection, with a significant impact on reservoir storage capacity.25 Therefore, a cushion gas with higher compressibility is required to allow the storage of large amounts of working gas with minimum pressure increase in the reservoir.26 Further, differences in viscosity and Received: February 20, 2019 Revised: June 7, 2019 Published: June 24, 2019 6527

DOI: 10.1021/acs.energyfuels.9b00518 Energy Fuels 2019, 33, 6527−6541

Article

Energy & Fuels

Figure 1. CH4 density (a), CO2 density (b), CH4 dynamic viscosity (c), and CO2 dynamic viscosity (d) changes with temperature and pressure.32,33 The black dots are the data points used for interpolation.

pressure evolution and gas mixing. Furthermore, the hydromechanical model is employed to study scenarios on the influence of the dip of geological structures as well as variations in reservoir thickness and temperatures on the mixing process. Finally, an assessment of reservoir integrity is undertaken by investigating changes in the reservoir’s elastic stability during the gas storage process. From that, a correlation between injection time and the reservoir’s tensile strength is obtained.

density between the working and cushion gases are required to mitigate gas mixing processes in the reservoir. Based on numerical simulations, Oldenburg26 discussed the physical properties of CO2 and CH4, and demonstrated that CO2 is likely to represent a suitable cushion gas due to its high effective compressibility. However, the authors investigated gas migration and mixing while neglecting the impact of accompanying temperature and geomechanical effects. Kühn et al.27,28 integrated subsurface CO2 and CH4 storage with the PtG approach, and studied the separate storage of CH4 and CO2 at depths, depending on general pressure and temperature gradients. Their focus was on the assessment of the economic feasibility of the integrated PtG and subsurface gas storage, but not the assessment of storing both gases in the same reservoir, or the extent of their mixing region. Niu and Tan29 studied the migration of gas mixtures using CO2 as the cushion gas in a subsurface gas storage reservoir, but neglected the dynamic changes of density and viscosity of CO2, CH4, and their mixtures at different pressure and temperature conditions. The purity of the working gas during gas withdrawal is particularly important and directly affects the energy efficiency of a PtG system. To study the migration and mixing of the working and cushion gases in the reservoir during the entire lifecycle of subsurface gas storage, a numerical multicomponent flow model considering mechanical processes and dynamic changes in fluid density and viscosity in the storage reservoir is presented. In this context, a multicomponent flow model is elaborated and successfully validated against published data26 based on the software package COMSOL Multiphysics.30 Since gas compressibility plays an important role in storage efficiency as demonstrated by Oldenburg,26 the reservoir compressibility as well as related porosity and permeability changes are expected to be important as well. Therefore, a hydromechanically coupled model is developed based on the validated multicomponent flow model to investigate the mechanical impacts on reservoir

2. MATERIALS AND METHODS 2.1. Physical Properties of CH4, CO2, and Their Mixtures. The density and dynamic viscosity of working and cushion gases affect the extent of the mixing region between the two gases in the reservoir. As shown in Figure 1a,c, the density and dynamic viscosity of CH4 increase steadily from 0 to 220 kg/m3 and from 1.2 × 10−5 to 2.6 × 10−5 Pa·s with increasing temperature and pressure, respectively. The steadiness of these changes results from the pressure and temperature conditions which are well above the critical point of CH4 (191.05 K, 4.54 MPa).31 Therefore, CH4 phase changes do not have to be considered in the multicomponent flow model. The situation for CO2 is different, since phase changes occur at temperatures above the critical temperature of 304.2 K and pressures higher than the critical pressure of 7.38 MPa. This occurrence induces substantial changes in density and dynamic viscosity, which must be represented in the numerical multicomponent flow model. As shown in Figure 1b,d, the density of CO2 changes in a range between 10 and 1000 kg/m3, and the dynamic viscosity changes between 0.2 × 10−4 and 1.4 × 10−4 Pa·s for the given temperature and pressure conditions. It is worth noting that CO2 density increases by a factor of 6 when the temperature increases from 300 to 320 K and the pressure from increases from 6 to 8 MPa. This illustrates the high CO2 compressibility in the vicinity of the critical point. For the same conditions, the density of CH4 increases only by a factor of 2. The dynamic viscosity of CO2 is approximately 10 times higher than that of CH4 at pressure and temperature conditions near the critical point of CO2. Although a significant difference in the densities and dynamic viscosities of the working and cushion gases should limit the gas 6528

DOI: 10.1021/acs.energyfuels.9b00518 Energy Fuels 2019, 33, 6527−6541

Article

Energy & Fuels

Figure 2. Density and dynamic viscosity of CO2−CH4 mixtures at temperatures of 303.15 K (a). Effect of temperature at 303.15, 313.15, and 323.15 K (b) on density and dynamic viscosity of gases.35,36

Figure 3. Sketch of a two-dimensional axisymmetric reservoir model and boundary conditions. mixing,34 the interface between CO2 and CH4 will tend to be unstable due to the high mobility ratio. As shown in Figure 2a, the physical properties of mixtures strongly correlate with the mole fraction of CH4. The physical properties of CO2 change substantially even if minor amounts of CH4 are added to the gas phase, whereby the magnitude of the change gradually decreases with increasing CH4 mole fraction. In addition, the physical properties of mixtures at temperatures of 303.15, 313.15, and 323.15 K are studied. To demonstrate the impacts of temperature on the physical properties, density and dynamic viscosity changes for CO2 and CH4 are provided

as a function of the CH4 mole fraction in the CO2 gas phase in Figure 2b. Here, only the pure CO2 and CH4 situations are shown because the effect of mole fraction is similar to that of Figure 2a. With the increasing amount of CH4 mixed with CO2, the density and viscosity of the mixed gases will decrease. Figure 2b illustrates the effect of temperature. It is apparent that temperature has little effect on the physical properties of CH4, but a notable impact on those of CO2. Further, Figure 2 shows that the effect of temperature gradually decreases with increasing CH4 mole fraction. For a temperature of 303.15 K in Figure 2a, the major changes in density and dynamic 6529

DOI: 10.1021/acs.energyfuels.9b00518 Energy Fuels 2019, 33, 6527−6541

Article

Energy & Fuels

where Ji is the diffusive mass flux of a single component in a multicomponent gas phase mixture (kg/(m2·s)); D is the diffusion coefficient (m2/s); ci is the concentration of a single component in the multicomponent gas phase mixture (kg/m3), where ci = χiρi; ρi·∇χi is

viscosity occur near the critical pressure, since the temperature is subcritical. However, as shown in Figure 2b, the changes in density and dynamic viscosity are smaller at higher temperatures. 2.2. Model Validation against Published Data. Since the commercial software COMSOL Multiphysics can be used to solve partial differential equations and obtain the solutions for coupled multiphysics problems,37 a multicomponent flow model is implemented using that software, and subsequently validated against Oldenburg’s modeling results.26 Four general simplifications are introduced into the mathematical model: (1) a single injection and production well is used; (2) the presence of water in the reservoir is neglected, since most of the groundwater is assumed to be displaced by the CO2 cushion gas; (3) the miscibility between CH4 and CO2 is taken into account based on the diffusion equation; and (4) the absorption effects of porous media are neglected. The model consists of a two-dimensional reservoir with a thickness of 22 m and a lateral extent of 1000 m (Figure 3). The vertical axis of symmetry at the injection and production well is used to reduce the model size to a half-domain model. The domain is discretized into 4400 uniform grid cells (200 × 22 grids in the X and Y directions, respectively). The injection well is located at the left model boundary at Y = −3 m (Figure 3). Neumann “no-flow” conditions are set at all model boundaries, and the initial pressure in the reservoir is constant at 6 MPa. The reservoir is initially fully saturated with CO2 at a temperature of 313.15 K. The injection period is 180 days. The parameters are consistent with Oldenburg’s model.26 To quantify the extent of the mixing zone over the entire reservoir, the reservoir midline is used as the reference (Y = −11 m). The mole fraction of CH4 in a range of 0.9−0.1 in the gas phase is regarded as the mixing region in the present study to allow for a comparison with Oldenburg’s modeling results. The distances dmin and dmax represent the nearest and farthest extents of the mixing region at the depth of the reference line calculated from the left boundary, respectively. The mixing degree of the two gases is defined by eq 1. η=

Vmix V

molecular diffusion; χi

∂ρi ∂T

∇T is

(5)

where ci = χiρi, and χi is the mole fraction of a gas component in the gas phase. In addition, the closure relation χCO2 + χCH4 = 1 is used. The two-phase Darcy flow module of the porous media and subsurface flow module in COMSOL Multiphysics are used with some modifications to quantify the mixing of CO2 and CH4 in the subsurface gas storage reservoir. First, the Darcy law equation in the two-phase Darcy flow module in COMSOL Multiphysics is modified to consider gravitational effects as shown in eq 3. Second, the density and dynamic viscosity of both gases are considered as functions of pressure, temperature, and CH4 mole fraction, obtained by linear interpolation of all data (Figures 1 and 2). All parameters used in the model are listed in Table 1.

Table 1. Properties of Gases and Reservoir Rocks

(2)

where ρ is density of the fluid mixture (kg/m3), which is a function of pressure (p), mole fractions of CH4 (xCH4) and CO2 (xCO2), and temperature (T) as shown in Figure 2; ϕ is the effective reservoir porosity; v is the Darcy velocity (m/s) of the fluid, which is derived from Darcy’s law (eq 3).

property

value

reservoir size initial porosity, ϕ0 initial permeability (isotropic), k0 diffusion coefficient, Dc temperature, T initial pressure, p0 CH4 injection rate, m0 elastic modulus, E Poisson’s ratio, υ rock density, ρr Biot−Willis coefficient, αp porosity, ϕ permeability, k density of CO2, ρCO2

22 × 1000 0.3 1 × 10−12 1 × 10−6 313.15 6 1.8375 × 10−2 3 0.3 1900 1 eq 6 eq 7 calcd in section 2.1

units

m2 kg/m3

dynamic viscosity of CO2, μCO2

calcd in section 2.1

Pa·s

density of CH4, ρCH4

calcd in section 2.1

kg/m3

dynamic viscosity of CH4, μCH4

calcd in section 2.1

Pa·s

density of gas mixtures, ρmix dynamic viscosity of mixtures, μmix

calcd in section 2.1 calcd in section 2.1

kg/m3 Pa·s

m m2 m2·s−1 K MPa kg·s−1 GPa kg/m3

2.3. Hydromechanically Coupled Model. The porosity and permeability of the reservoir are subject to changes with increasing pore pressure during the injection process, which also affects fluid migration. A specific model is built based on the multicomponent flow model to study the hydromechanical interaction resulting from gas storage, with a special focus on the mixing process. Porous reservoir rock is considered as a linearly elastic material. Based on the generalized Hooke’s law, the control equation of the stress field is expressed by displacement and pore pressure p, which can be written as eq 6.

(3)

where k is the reservoir permeability (m ); μ is the viscosity of the fluid mixture (Pa·s), which is a function of pressure (p), mole fractions of CH4 (xCH4) and CO2 (xCO2), and temperature (T) as shown in Figure 2. p is the pore pressure (Pa), and g is the gravitational acceleration (9.81 m/s2). According to Fick’s law, the diffusive flux of a fluid component is proportional to its concentration.38 2

∂ρ ji Ji = − D∇ci = − D(ρi ·∇χi + χi ·∇ρi ) = − Djjjρi ·∇χi + χi i j ∂p k yz ∂ρ ∇p + χi i ∇T zzz z ∂T {

∇p is barodiffusion; χi

∂(ϕci) + ∇(ci v) = ∇(Dc∇ci) ∂t

(1)

k v = − (∇p + ρg ) μ

∂p

thermodiffusion. The diffusion equation (eq 5) can be obtained by combining eq 4 with eq 2.

where η is the mixing degree, Vmix is the volume of the mixing area, and V is the volume of the reservoir. The spatial distribution and distance of the mixing region extent in the reservoir are studied after 30, 90, and 180 days of CH4 injection. The continuity equation (eq 2) for a single phase with two gas components is based on the mass conservation equations and is used to estimate the flow field in the model.

∂(ρϕ) + ∇(ρv) = 0 ∂t

∂ρi

Gui , jj + (G + λ)uj , ij + αppi + Fi = 0 where G =

E 2(1 + υ)

(6)

is the shear modulus (MPa); ui is the displacement

(m) component in the ith direction (i = x, y); λ =

Eυ (1 + υ)(1 − 2υ)

is the

Lamé constant (MPa); ui,jj is the shear stress component in the jth direction (j = y, x); uj,ij is the normal stress component in the jth direction (j = y, x); αp is the Biot−Willis coefficient; pi is the pore

(4) 6530

DOI: 10.1021/acs.energyfuels.9b00518 Energy Fuels 2019, 33, 6527−6541

Energy & Fuels

ÄÅ ÉÑ3 k 0 ÅÅÅÅ εv Δp ÑÑÑÑ + (1 − ϕ0) k= Å1 + Ñ ϕ0 ϕ0K s ÑÑÑÑÖ 1 + εv ÅÅÅÅÇ

pressure (MPa); Fi = ρrg is the volumetric force (kg/(m·s)2) in the ith direction (i = x, y); ρr is the density of the reservoir (kg/m3); E is the elastic modulus of the reservoir (MPa); υ is Poisson’s ratio. The hydromechanical modeling process is shown in Figure 4. The pore pressure in the flow field will affect stresses and strains in the

Article

(8)

where εv is the volumetric strain and εv = εx + εy + εz; εij = (1/2)(ui,j + uj,i); ϕ0 is the initial reservoir porosity; Δp is the increment in pore pressure (MPa); Ks is the solid grain bulk modulus (MPa); k0 is the initial permeability (m2). The Solid Mechanics module in COMSOL Multiphysics is coupled with the multicomponent flow model. Porosity and permeability are defined by the dynamic evolution model shown in eqs 7 and 8 and written as variables in the component definitions part in COMSOL Multiphysics. The applied physical properties of both fluids and the reservoir are listed in Table 1. 2.4. Impacts of Reservoir Thickness, Formation Dip, and Reservoir Temperature. The thicknesses of geological reservoirs may vary from a few to hundreds of meters. To quantify the influence of the reservoir thickness on the spatial extent of the mixing region, reservoir thicknesses of 22, 50, and 100 m are studied using the hydromechanical model. The mixing degree is obtained by eq 1, and the angle between the miscible interface and horizontal plane is used to describe the position of the mixing region. All other parameters assigned to the model are listed in Table 1. Further, geological reservoirs may have different structural geological features, including varying dip angles, which likely have an impact on subsurface gas storage when CO2 is used as the cushion gas. To quantify these impacts, reservoirs with dip angles of ±5 and ±10° have been employed in the hydromechanically coupled model to represent syncline and anticline structures, respectively. The mixing degree and location of the mixing region are studied in these scenarios with all other parameters listed in Table 1.

Figure 4. Hydromechanical modeling framework.

reservoir due to Therzaghi’s effective stress principle. Variations in stresses and strains will induce displacements in the reservoir and affect its porosity and permeability. The hydromechanically coupled model is obtained by alternating the mechanical and hydraulic variables in COMSOL Multiphysics. The mechanical calculations are affected by the presence of a fluid pressure in the pores, resulting in deformation of the porous media due to effective stresses. The flow calculations are affected by poroelastic porosity and permeability changes during the gas injection. Those two reservoir parameters are employed to describe fluid flow as expressed in eqs 7 and 8.39 ÅÄ ÑÉ Δp ÑÑÑ 1 ÅÅÅÅ ϕ= ÅÅϕ0 + εv + (1 − ϕ0) ÑÑÑ 1 + εv ÅÅÇ Ks ÑÑÖ (7)

Figure 5. Pressure evolution at injection point in Oldenburg’s model,26 multicomponent flow model and hydromechanically coupled model with CO2 (a) and CH4 (b) cushion gases. 6531

DOI: 10.1021/acs.energyfuels.9b00518 Energy Fuels 2019, 33, 6527−6541

Article

Energy & Fuels

Figure 6. Distance of mixing region from the injection well after 30, 90, and 180 days of CH4 injection. Reservoir temperature affects the densities and dynamic viscosities of CH4, CO2, and their mixtures. To quantify the impact of temperature on the spatial extent of the mixing region, temperatures of 303.15, 313.15, and 323.15 K are used, while all other parameters remain unchanged. The mixing degree and pressure evolution in the reservoir are quantified at these temperatures. 2.5. Reservoir Stability. Because the modeling process is based on the linear elastic theory, the results are only valid as long as the reservoir behaves elastically in the absence of shear and tensile failure. Therefore, the mechanical stability of the reservoir is monitored during the injection process. The von Mises failure criterion is used here to evaluate the risk of reservoir failure (eq 9). g(σ ) =

(σp1 − σp2)2 + (σp2 − σp3)2 + (σp1 − σp3)2 2σts 2

time-dependent pressure evolution in Oldenburg’s model and our multicomponent flow model is 0.9970, with a mean difference of 0.16 MPa. Hereby, the average difference in relation to Oldenburg’s model is 2.01%. In Figure 5b, the coefficient of determination between the pressure data in Oldenburg’s model and our multicomponent flow model is 0.9996, with a mean difference of 0.09 MPa, whereby the average difference in relation to Oldenburg’s model is 0.9%. The spatial position and the extent of mixing region are studied and compared in the following. As shown in Figure 6, the migration distance changes substantially after 30, 90, and 180 days of CH4 injection. In the multicomponent flow model, dmin and dmax are 116.9 m and 234.2 m after 30 days of injection, 342.7 m and 568.3 m after 90 days of injection, and 530.3 m and 786.8 m after 180 days of injection, respectively. The average migration distance is 175.5, 455.5, and 658.5 m after 30, 90, and 180 days of injection, respectively. The location of the mixing region in the multicomponent flow model almost coincides with that provided by Oldenburg’s simulation results.26 Based on numerical simulation and data postprocessing, the degrees of mixing after 30, 90, and 180 days are shown in Table 2.

−1≥0

(9) where g(σ) is the safety factor, whereby g(σ) ≥ 0 represents reservoir failure; σp1 is the maximum principal stress; σp2 is the intermediate principal stress; σp3 is the minimum principle stress; σts is the tensile strength. Based on eq 9, the reservoir equivalent stress intensity σ̅ can be obtained from the hydromechanical model, where σ̅ =

1 [(σp1 − σp2)2 + (σp2 − σp3)2 + (σp1 − σp3)2 ] 2

In other words, the mechanical stability of the reservoir is maintained if the tensile strength remains below the equivalent stress intensity. All the parameters used in the respective simulation runs are listed in Table 1. The equivalent stress intensity of the reservoir is studied for injection rates of 1.8375 × 10−2 and 4.00 × 10−2 kg/s/m. At these injection rates, reservoir integrity is not compromised if the tensile strength of reservoir is assumed with 5 MPa.

Table 2. Comparison of Values of Mixing Degree after 30, 90, and 180 days of Injection in the Investigated Scenarios η (%)

3. RESULTS 3.1. Multicomponent Flow Model Validation. The pressure evolution in the scenarios using CO2 and CH4 as cushion gases are shown in Figure 5, parts a and b, respectively. The pore pressure increases significantly if methane is used as cushion gas, while a slower increase is observed with a CO2based cushion gas. Figure 5 outlines the pressure evolution produced by Oldenburg’s model26 as well as the multicomponent flow and hydromechanical model developed in the present study. In Figure 5a, the coefficient of determination (R2) between the

time (days)

Oldenburg’s model

multicomponent flow model

hydromechanically coupled model

30 90 180

7.89 13.57 15.65

6.69 14.21 18.47

5.21 12.89 16.73

Since two-dimensional models are considered in the present study, the mixing area ratio η is an effective way to quantify the mixing degree between both gases. Table 2 shows that, with an increase in injection time, the mixing degree increases significantly. By comparing the results of Oldenburg’s model26 and our multicomponent flow model, the deviation is less than 3%, indicating that our modeling approach is valid and sufficiently accurate to represent the relevant processes. 6532

DOI: 10.1021/acs.energyfuels.9b00518 Energy Fuels 2019, 33, 6527−6541

Article

Energy & Fuels 3.2. Hydromechanically Coupled Model. In the hydromechanical model, pore pressure increases to 11.3 MPa when CH4 is employed as cushion gas, while it peaks at 9.37 MPa if CO2 is used for that purpose after 180 days of injection (Figure 5). It is worth noting that pressure evolution in the hydromechanical model is lower than that in the multicomponent flow model. The maximum differences of pore pressure between our hydromechanical and Oldenburg’s model are 2.08 MPa for the CO2 cushion gas and 2.41 MPa for the CH4 cushion gas after 180 days of injection. As shown in Figure 6, dmin and dmax calculated by the hydromechanical model are 127.5 m and 193.9 m after 30 days of injection, 311.4 m and 458.2 m after 90 days of injection, and 487.0 m and 751.8 m after 180 days of injection, respectively. The average distances are 160.7, 384.8, and 619.4 m after 30, 90, and 180 days of injection, respectively. Compared to the multicomponent flow model, the mixing region migration distance in the hydromechanical model is smaller. The values of the mixing degree shown in Table 2 are 1.48, 1.32, and 1.74% smaller than those produced by the multicomponent flow model after 30, 90, and 180 days of injection. The pressure evolution at the points (100,0), (500,0), and (900,0) along the reservoir is studied. The result shows that, after the same time of injection, a point far from the injection point has a smaller pressure than one closer to the injection point. Because the flow boundary conditions in the reservoir are closed, the pressure difference is quite small (less than 0.5 MPa). The density and dynamic viscosity of the mixed gases are shown in Figure 7b,c. After 8, 68, and 165 days of injection, the points at (100,0), (500,0), and (900,0), respectively, begin to mix. The porosity increased by 0.002, and the permeability increased by 0.2 × 10−13 m2. 3.3. Impact of Reservoir Thickness. The distribution of the mixing region is again calculated for 30, 90, and 180 days of CH4 injection based on the hydromechanical simulation results shown in Figure 8a. The legend shown on the right side gives the mole fraction of CH4, with values ranging from 0 to 1. The white arrows represent the velocity magnitude and direction of fluid flow in the reservoir. The velocity of fluid flow in the X direction after 30 days of injection is 5 × 10−5 m/s, and gradually decreases with increased injection duration. The inclination angle of the miscible interface is approximately 26° after 180 days of injection. The mixing region distribution after 30, 90, and 180 days of CH4 injection without considering the gravity effects is shown in Figure 8b. The miscible interface is nearly perpendicular to the reservoir during the injection process. In this case, the degree of mixing is obviously smaller than a case that considers gravity effects. However, the failure to consider gravity yields an overly simple, ideal model that neglects the effect of the huge density difference between the two kinds of gases. There will always be a horizonal flow in the reservoir. The results may not provide a good reference for practical application. Therefore, gravity should be considered when studying the impact of reservoir thickness. With increasing reservoir thickness, as shown in Figure 8c,d, the degree of mixing is reduced. At a reservoir thickness of 50 m, the mixing degrees are 2.28, 6.59, and 13.29% after 30, 90, and 180 days of injection, respectively. When the reservoir thickness is assumed to be 100 m, the mixing degree is smaller, with 0.90, 2.58, and 5.57% after 30, 90, and 180 days of

Figure 7. Pressure evolution (a), density (b), and dynamic viscosity (c) of gases at (100,0), (500,0), and (900,0) along the reservoir; reservoir porosity (d) and permeability (e) change during the injection process in hydromechanically coupled model.

injection, respectively. The inclination angle of the miscible interface is smaller and tends in the horizontal direction. When the reservoir thickness is set to 50 m, the inclination angle of the miscible interface is approximately 16° after 180 days of 6533

DOI: 10.1021/acs.energyfuels.9b00518 Energy Fuels 2019, 33, 6527−6541

Article

Energy & Fuels

Figure 8. (a) Distribution of mixing region in a 22 m thick reservoir after injection periods of 30, 90, and 180 days calculated with the hydromechanically coupled model. (b) Distribution of mixing region in a 22 m thick reservoir after injection periods of 30, 90, and 180 days calculated with the hydromechanically coupled model without considering gravity effects. (c) Distribution of mixing region in 50 m thick reservoir after 30, 90, and 180 days of injection, calculated by the hydromechanical model. (d) Distribution of mixing region in a 100 m thick reservoir after 30, 90, and 180 days of injection calculated by the hydromechanical model.

injection. This angle is reduced to approximately 11° for a reservoir thickness of 100 m. For the present axisymmetric model, the distribution region of CH4 gradually becomes an inverted cone with increasing reservoir thickness, which is particularly evident in Figure 8d. 3.4. Effects of Structural Geological Features. To quantify the impact of structural geological features, i.e., the dip angle, on the mixing region, syncline and anticline structures with different dip angles are studied. The mixing region distribution in the reservoir at dip angles of 5 and 10° are shown in Figure 9. For viewing purposes, the Y-axis is

scaled down to 1/3 of the original length. This model is used to represent a syncline structure. Figure 9 shows that the mixing region in the syncline structure is much larger than that in the horizontal reservoir presented in Figure 8a. In a syncline reservoir with a 5° dip, the degree of mixing is 5.46, 15.79, and 30.13% after 30, 90, and 180 days of injection, respectively. For a 10° dip, the degree of mixing is 5.71, 20.17, and 40.81% after 30, 90, and 180 days of injection, respectively. The mixing region distribution increases with increasing dip angles (cf. Figure 9). After 30 days of injection, the mixing region distribution ranges from 0 m to approximately 280 m in the reservoir with a dip angle of 5°, 6534

DOI: 10.1021/acs.energyfuels.9b00518 Energy Fuels 2019, 33, 6527−6541

Article

Energy & Fuels

Figure 9. (a) Distribution of mixing region in a syncline reservoir with 5° dip after 30, 90, and 180 days of injection. (b) Distribution of mixing region in a syncline reservoir with 10° dip after 30, 90, and 180 days of injection.

but from 0 m to approximately 320 m in the reservoir with a dip angle of 10°. After 90 days of injection, the difference in the distribution of the mixing region is substantially larger, ranging from approximately 140 to 650 m in the reservoir with a dip angle of 5°, but from approximately 70 to 720 m in the reservoir with a dip angle of 10°. After 180 days of injection, a notable mixture area develops at the right boundary of the reservoir for a dip angle of 10°. The mixing region distribution in the anticline reservoir with dip angles of 5 and 10° are shown in Figure 10. Compared with the results presented in Figures 8a and 9a,b, the mixing region is the smallest in the anticline structure at the same simulation time. In the 5° dip syncline reservoir, the degree of mixing is 3.86, 8.56, and 12.48% after 30, 90, and 180 days of injection, respectively. In the 10° dip syncline reservoir, the degree of mixing is 3.27, 4.23, and 5.74% after 30, 90, and 180 days of injection, respectively. Further, the distribution of the mixing region decreases with increasing dip angles. After 30

days of injection, the distribution of the mixing region ranges from approximately 80 to 240 m in the reservoir with a dip angle of 5°, while it only ranges from approximately 100 to 210 m in the reservoir with a dip angle of 10°. After 90 and 180 days of injection, the difference in the mixing region distribution in reservoirs with dip angles of 5 and 10° is substantially larger. The extent of the mixing region at a reservoir dip angle of 10° is approximately 80 m smaller than for that with a dip angle of 5° after 90 days of injection. After 180 days of injection, this difference amounts to approximately 120 m. 3.5. Influence of Reservoir Temperature. Figure 11 shows that the mixing region tends to decrease with increasing temperature. To clearly show the differences, the calculated mixing degree η is presented in Table 3. After 30 days of injection, the degree of mixing is 4.98% at 323.15 K, which is smaller compared to those at 303.15 and 6535

DOI: 10.1021/acs.energyfuels.9b00518 Energy Fuels 2019, 33, 6527−6541

Article

Energy & Fuels

Figure 10. (a) Distribution of mixing region in an anticline reservoir with dip angle of 5° after 30, 90, and 180 days of injection. (b) Distribution of mixing region in an anticline reservoir with 10° dip after 30, 90, and 180 days of injection.

Figure 11. Distribution of mixing region at 303.15, 313.15, and 323.15 K after 30, 90, and 180 days of injection.

313.15 K. This trend prevails after 90 and 180 days of injection. In the hydromechanical model, the pore pressure evolution at different temperatures is also studied. As shown in Figure

12, the reservoir pore pressure is lower at 303.15 K than that at 313.15 and 323.15 K. The pore pressure increases to 8.5, 9.4, and 10.0 MPa after 180 days of injection at 303.15, 313.15 and 323.15 K, respectively. 6536

DOI: 10.1021/acs.energyfuels.9b00518 Energy Fuels 2019, 33, 6527−6541

Article

Energy & Fuels

We obtained those by linearly interpolating data from Kunz and Wagner35 and NIST36 for our multicomponent flow model, while Oldenburg’s model uses the Peng−Robinson equation of state implemented in TOUGH2 EOS7C.40 The increase in pore pressure is lower when CO2 is used as the cushion gas compared to the CH4 case. Consequently, more CH4 can be stored in the reservoir at the same pressure conditions if CO2 is used as the cushion gas. It is worth noting that only a small amount of CH4 mixed with the CO2 cushion gas results in an immediate decrease in the mixture’s density. For that reason, the weight-average model incorporated in COMSOL Multiphysics cannot be applied to represent the physical properties of the gas mixtures. A hydromechanically coupled model was implemented based on the previously introduced multicomponent flow model, considering the interaction between hydraulic and mechanical effects, dynamic porosity, and dynamic permeability changes. As shown in Figures 5 and 6, the pressure evolution and distance of mixing region migration are both lower in our hydromechanical model than in Oldenburg’s model or the multicomponent flow model at identical injection times. Because of the effective stress principle, reservoir deformation lowers reservoir pore pressure, which is beneficial for increasing the storage capacity of the working gas. At decreasing pore pressures, the mixing region migration distance is also decreasing. In addition, the densities of CO2, CH4, and their mixtures are reduced due to the decreased pore pressure in the reservoir, so that the concentration of the single components in the multicomponent gas mixture (ci = χiρi) is also reduced. Based on eq 4, the diffusive mass flux of a single component in multicomponent gas mixture (Ji) is consequently decreased as well. Therefore, the degree of mixing in the hydromechanical model is notably lower than that in the multicomponent flow model. However, the degree of mixing in the hydromechanical model is higher than that in Oldenburg’s model after 180 days of injection. This is because reservoir porosity and permeability are not equal to their initial hydromechanical model values at that time. Reservoir deformation increases both parameters, and then increases the fluid flow velocity based on Darcy’s law, which gradually increases their miscibility. Geological reservoirs with different thicknesses have a significant impact on the distribution of both gases in the reservoir following the CH4 injection. In this context, applying the hydromechanical model to assess the behavior of the

Table 3. Mixing Degree (%) at 303.15, 313.15, and 323.15 K after 30, 90, and 180 days of Injection time (days) temp (K)

30

90

180

303.15 313.15 323.15

6.25 5.21 4.98

15.52 12.89 11.48

19.81 16.73 15.73

3.6. Mechanical Stability of the Reservoir. Since all calculations are based on linear elasticity, the simulation results become inaccurate when the reservoir is stressed beyond the elastic limits due to the occurrence of shear and/or tensile failure. Therefore, the development of the reservoir’s mechanical stability during the injection process is taken into account in the following. The hydromechanical model in the present study can describe the correlation between days of injection without any occurrence of failure and the tensile strength of the reservoir. In other words, the mechanical stability of reservoir can be maintained if the tensile strength of the reservoir remains below the equivalent reservoir stress intensity. The equivalent stress intensity of the reservoir varies from 3.43 to 5.96 MPa at an injection rate of 1.8375 × 10−2 kg/s for 180 days of injection (Figure 13). With an increased injection rate, the required equivalent stress intensity of the reservoir increases. If the tensile strength of the reservoir rocks is assumed to be 5 MPa, the first rock failure will occur after 68 days at an injection rate of 4.00 × 10−2 kg/s and after 113 days at an injection rate of 1.8375 × 10−2 kg/s. Figure 14 shows the correlation between an injection time without any occurrence of failure and the reservoir’s tensile strength at an injection rate of 4.00 × 10−2 kg/s. These findings show that a higher tensile strength also prolongs the injection time without an occurrence of failure.

4. DISCUSSION The mixing region between CH4 and CO2 for subsurface gas storage using CO2 as the cushion gas was studied using the presented multicomponent flow model, which is validated against Oldenburg’s model.26 The comparison shows a very good agreement between the models, with a maximum deviation of less than 3%. The slight difference between the models likely arises from the equations of state used to calculate the physical properties of gases and their mixtures.

Figure 12. Pressure evolution at injection point for 303.15, 313.15, and 323.15 K. 6537

DOI: 10.1021/acs.energyfuels.9b00518 Energy Fuels 2019, 33, 6527−6541

Article

Energy & Fuels

Figure 13. Equivalent stress intensity of reservoir at two different injection rates.

Figure 14. Correlation between injection time without failure and tensile strength at injection rate of 4.00 × 10−2 kg/s.

syncline structures. Larger dip angles increase gravitational effects, and thus increase the degree of mixing. At the same time of CH4 injection, anticline structures show an obvious advantage in terms of mixing region extent control. CO2 migrates downward, which avoids secondary mixing with CH4, and preserves the purity of CH4 in the flooding area. Therefore, compared to syncline and horizontal storage reservoirs, the lowest mixing region extent occurs in anticline structures, and the mixing region can become smaller with increasing dip angles of anticline structures. A large anticline dip is preferred for underground gas storage that uses CO2 as a cushion gas. At a temperature of 303.15 K, the densities of CO2, CH4, and their mixtures are higher than those at 313.15 and 323.15 K, close to the critical CO2 pressure (Figure 2). Hence, the concentration of a single component in the multicomponent gas mixture ci is increased (ci = χiρi) so that the diffusive mass flux Ji will increase if a constant diffusion coefficient is assumed, based on eq 4. Therefore, mixing is decreased at higher temperatures. In terms of the pressure evolution shown in Figure 12, the pore pressure is lowest at a temperature of 303.15 K. Hereby, the compressibility of CO2 directly affects the pore pressure. As shown in Figure 15, the compression factor varies from 0 to 1, and represents the effort required to compress the gaseous phase. When the compression factor is 1, the fluid is incompressible. CO2 shows a high compressibility close to its critical pressure, whereby it becomes less compressible with increasing temperature. The equivalent stress intensity of the reservoir obtained by using hydromechanical simulations has to be less than the

mixing region with different conditions is highly relevant for the efficient operation of a gas storage site. It has been demonstrated that CO2 is displaced and compressed at the reservoir boundaries with the continuous injection of CH4. Due to the relatively high CO2 compressibility close to its critical point, it represents an effective cushion gas. After 90 and 180 days of injection, the fluid flow velocity gradually decreases with continuous injection. This is due to the decrease in CO2 compressibility with increasing pore pressure. It is worth noting that the degree of mixing and inclination angle of the miscible interface decrease with increasing reservoir thickness. This can be explained by the available pore space in a thick reservoir, which provides sufficient volume for the gravitational separation of both fluids. The high differences in density and dynamic viscosity between CO2 and CH4 (Figure 2) are the key factors contributing to gas separation. A thick reservoir will transform the horizontal flow into upper and lower stratified flows due to the high density difference in the reservoir. Meanwhile, the fluid flow velocity at reservoir thicknesses of 50 and 100 m is significantly reduced due to the lower pore pressure increase. This also indicates that the reservoir still has excess capacity for further gas storage. Geological structures have a high impact on gas mixing processes in the reservoir. The degree of mixing is especially high in syncline reservoirs and becomes even higher at increasing dip angles. At the same temperature and pressure conditions, the density of CO2 is significantly higher than that of CH4. During the injection process, backward flow CO2 will occur due to gravity effects, so the mixing area is enlarged in 6538

DOI: 10.1021/acs.energyfuels.9b00518 Energy Fuels 2019, 33, 6527−6541

Article

Energy & Fuels

the kinetic theory, diffusion is more intensive for gases consisting of small molecules. For instance, the diffusion of hydrogen in air is 5 times higher than that of oxygen.38 In our model, the diffusion coefficient is considered for the CH4− CO2 mixing system, which is set as 1 × 10−6 m2 s−1, as same as the value from Oldenburg’s model. It is true that the diffusion coefficient we set in our model is convincing to represent the gases in the mixing area, rather than all the gases in the reservoir. The self-diffusion coefficients of CO2 and CH4 should be considered as well for further study.

5. CONCLUSIONS The main conclusions from the aforementioned investigation are summarized as follows: A multicomponent flow model was implemented and validated against Oldenburg’s model, with a deviation of less than 3%, indicating that our model complies with the purpose of the current study. A hydromechanical model was developed based on the multicomponent flow model to consider compressibility effects in the reservoir. In the hydromechanical model, the pore pressure, mixing region migration distance, and degree of mixing are lower than those in the multicomponent flow model, which indicates beneficial effects in terms of the working gas storage capacity of the reservoir with reduced mixing effects. Reservoirs with large thicknesses tend to have lower pore pressure, and thus smaller mixing region migration, which provides sufficient time and space for gas separation due to the high difference in density and dynamic viscosity between CO2 and CH4. Thus, a reservoir with a higher thickness is better suited to reduce the extent of the mixing region. The influence of structural geological features on the distribution of the mixing region shows that it has the lowest extent in anticline structures, compared to syncline and horizontal storage reservoirs. In addition, the mixing region distribution decreases with increasing dip angles in anticline structures. Therefore, large anticline dips are preferred for underground gas storage that uses CO2 as a cushion gas. The degree of mixing is the lowest at a temperature of 323.15 K in our study. However, it decreases with increasing temperature in our hydromechanical simulations, since the reservoir pore pressure increases with increasing temperature. Reservoir stability is strongly related to the reservoir’s tensile strength and is maintained if the tensile strength remains below the reservoir equivalent stress intensity. A safe injection rate can be obtained from the stability analysis presented in this study.

Figure 15. Compression factor of CO2 at 303.15, 313.15, and 323.15 K.

tensile strength of the reservoir to maintain a safe injection process. With a certain CH4 injection rate, a relationship between the injection time without any occurrence of failure and the reservoir’s tensile strength is given to enable the applicability of the results to gas storage operation. However, groundwater flow and preexisting fractures can have a substantial impact on the gas flow in a gas storage reservoir, factors which will be considered in future research activities. Furthermore, feasible gas injection rates, injection pressures, and injection well locations will be the assessed in further studies. Compared with Oldenburg’s model and our multicomponent flow model, hydromechanical coupling is required to obtain more accurate results on the impact of gas mixing in subsurface storage reservoirs. Due to the strongly coupled interactions between mechanical effects, reservoir pressure evolution and pore pressures in the hydromechanically coupled model are lower for gas mixture systems. The distance of the mixing region from the injection well and the degree of mixing are smaller than those calculated by Oldenburg’s and our multicomponent flow models. These results are especially relevant for subsurface gas storage site operators, and must be considered to determine the economics of gas storage. However, there are still some shortcomings in the present study to be addressed in future studies: we have employed an isothermal modeling approach. Further, single-phase multicomponent flow has been considered assuming that groundwater has been already displaced from the reservoir by the cushion gas. There still may be further impacts on the pressure development in the reservoir in case an aquifer model is taken into account instead of no-flow boundaries. The effect of gas diffusion on the mixing after gas injection, especially during long-time storage, is not studied. It is very likely that the mixing area will continue to expand due to gas diffusion effects. An optimized recovery method is necessary to reduce the mixing during the withdrawal period. In addition, our reservoir mechanical stability analysis is limited to tensile failure following the von Mises criterion, due to the linear elasticity assumption. It is worth noting that all the results are based on the same diffusion coefficients of CH4 and CO2. The diffusion coefficient can be divided into the self-diffusion coefficient of a certain kind of gas41 and the diffusion coefficient of gases in a mixing system (such as CH4−CO2 mixing system). Based on



AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Tel.: +86-27-87198126. Fax: +8627-87198967. ORCID

Qi Li: 0000-0003-0679-4385 Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS This work is mainly supported by IRSM-GFZ Subsurface Utilization of Captured Carbon and Energy Storage System 6539

DOI: 10.1021/acs.energyfuels.9b00518 Energy Fuels 2019, 33, 6527−6541

Article

Energy & Fuels

(2) Celia, M. A. Geological storage of captured carbon dioxide as a large scale carbon mitigation option. Water Resour. Res. 2017, 53 (5), 3527−3533. (3) Monastersky, R. Global carbon dioxide levels near worrisome milestone. Nature 2013, 497 (7447), 13−4. (4) Vishal, V.; Ranjith, P.; Singh, T. CO2 permeability of Indian bituminous coals: Implications for carbon sequestration. Int. J. Coal Geol. 2013, 105, 36−47. (5) Peters, G. P.; Andrew, R. M.; Boden, T.; Canadell, J. G.; Ciais, P.; Le Quéré, C. L.; Marland, G.; Raupach, M. R.; Wilson, C. The challenge to keep global warming below 2 °C. Nature Climate Change 2013, 3 (1), 4−6. (6) Zhang, Z.; Huisingh, D. Carbon dioxide storage schemes: Technology, assessment and deployment. J. Cleaner Prod. 2017, 142, 1055−1064. (7) Meinshausen, M.; Meinshausen, N.; Hare, W.; Raper, S. C.; Frieler, K.; Knutti, R.; Frame, D. J.; Allen, M. R. Greenhouse-gas emission targets for limiting global warming to 2 °C. Nature 2009, 458 (7242), 1158−1162. (8) Liu, H.; He, Q.; Borgia, A.; Pan, L.; Oldenburg, C. M. Thermodynamic analysis of a compressed carbon dioxide energy storage system using two saline aquifers at different depths as storage reservoirs. Energy Convers. Manage. 2016, 127, 149−159. (9) Chen, D.; Pan, Z.; Liu, J.; Connell, L. D. Modeling and simulation of moisture effect on gas storage and transport in coal seams. Energy Fuels 2012, 26 (3), 1695−1706. (10) Chandler, H. Empowering Variable Renewables Options for Flexible Electricity Systems; International Energy Agency: Paris, France, 2008. (11) Niu, S.; Wang, Y.; Qu, W.; Qiang, W. Wind Power Industry in West China: Difficulties and SolutionsReflections on “Abandonment of Wind Power”. Sci. Technol. Rev. 2017, 10, 11−12. (in Chinese with English abstract). (12) Götz, M.; Lefebvre, J.; Mörs, F.; McDaniel Koch, A.; Graf, F.; Bajohr, S.; Reimert, R.; Kolb, T. Renewable Power-to-Gas: A technological and economic review. Renewable Energy 2016, 85, 1371−1390. (13) Ma, J.; Li, Q.; Kühn, M.; Nakaten, N. Power-to-gas based subsurface energy storage: A review. Renewable Sustainable Energy Rev. 2018, 97, 478−496. (14) Kühn, M.; Li, Q.; Nakaten, N. C.; Kempka, T. Integrated subsurface gas storage of CO2 and CH4 offers capacity and state-ofthe-art technology for energy storage in China. Energy Procedia 2017, 125, 14−18. (15) Castellani, B.; Rinaldi, S.; Morini, E.; Nastasi, B.; Rossi, F. Flue gas treatment by power-to-gas integration for methane and ammonia synthesis−energy and environmental analysis. Energy Convers. Manage. 2018, 171, 626−634. (16) Nakaten, N.; Schlüter, R.; Azzam, R.; Kempka, T. Development of a techno-economic model for dynamic calculation of cost of electricity, energy demand and CO2 emissions of an integrated UCG−CCS process. Energy 2014, 66, 779−790. (17) Tek, M. R. Underground Storage of Natural Gas; Gulf Publishing Co.: 1989; pp 15−21. (18) Blanco, H.; Faaij, A. A review at the role of storage in energy systems with a focus on Power to Gas and long-term storage. Renewable Sustainable Energy Rev. 2018, 81, 1049−1086. (19) Vandewalle, J.; Bruninx, K.; D’haeseleer, W. Effects of largescale power to gas conversion on the power, gas and carbon sectors and their interactions. Energy Convers. Manage. 2015, 94, 28−39. (20) Misra, B. R.; Foh, S. E.; Shikari, Y. A.; Berry, R. M.; Labaune, F. The Use of Inert Base Gas in Underground Natural Gas Storage. SPE Gas Technology Symposium; Society of Petroleum Engineers: 1988. DOI: 10.2118/17741-MS. (21) Gümrah, N. K. Fevzi, A numerical simulation study on mixing of inert cushion gas with working gas in an underground gas storage reservoir. Energy Sources 2000, 22 (10), 869−879.

(Grant No. 1plus1-2018-01) and the China Scholarship Council (Grant No. 201804910875).



NOMENCLATURE ci = concentration of a single component in multicomponent gases (kg/m3) D = diffusion coefficient (m2/s) E = elastic modulus (Pa) Fi = volumetric force (kg/(m·s)2) g = gravitational acceleration (9.81 m/s2) G = shear modulus (Pa) Ji = diffusive mass flux (kg/(m2·s)) k = permeability (m2) k0 = initial permeability (m2) Ks = solid grain bulk modulus (Pa) p = pore pressure (Pa) R2 = coefficient of determination χi = mole fraction of a component gas T = temperature (K) Vmix = volume of mixing area (m3) V = volume of reservoir (m3) X = axial direction X Y = axial direction Y Z = axial direction Z

Abbreviations

CCS = carbon dioxide capture and storage CH4 = methane CO2 = carbon dioxide H2 = hydrogen NIST = National Institute of Standards and Technology N2 = nitrogen PtG = power-to-gas Greek Symbols

αp = Biot−Willis effective coefficient ρr = rock density (kg/m3) μ = fluid viscosity (Pa·s) ui = displacement (m) v = Darcy velocity (m/s) υ = Poisson’s ratio λ = Lamé constant (Pa) ρ = fluid density (kg/m3) ϕ = effective porosity εv = volumetric strain ϕ0 = initial porosity Δp = change in pore pressure (Pa) g(σ) = safety factor σp1 = maximum principal stress (Pa) σp2 = intermediate principal stress (Pa) σp3 = minimum principle stress (Pa) σts = tensile strength (Pa) σ̅ = equivalent stress intensity (Pa) xCH4 = mole fraction of CH4 η = area ratio, i.e., mixing degree Subscripts

i = direction, i = x, y, z mix = mixing



REFERENCES

(1) Mathias, S. A.; Hardisty, P. E.; Trudell, M. R.; Zimmerman, R. W. Screening and selection of sites for CO2 sequestration based on pressure buildup. Int. J. Greenhouse Gas Control 2009, 3 (5), 577−585. 6540

DOI: 10.1021/acs.energyfuels.9b00518 Energy Fuels 2019, 33, 6527−6541

Article

Energy & Fuels (22) Obro, H. Study of underground gas storage using nitrogen as cushion gas. Presented at the 17th World Gas Conference, Washington, DC, USA, June 5, 1988. (23) Jiao, W. L.; Wang, Z. S.; Juan-Juan, L. I. Simulation on use of inert cushion gas in underground gas storage reservoir in aquifer. J. Harbin Inst. Technol. 2008, 40 (12), 1949−1955. (in Chinese with English abstract). (24) Li, P.; Li, J. Numerical Simulation of Using Inert Gas as Cushion Gas in Aquifer Underground Gas Storage Reservoir. Gas Heat 2012, 32 (5), A35−A42. (in Chinese with English abstract). (25) Kim, J.; Choi, J.; Park, K. Comparison of nitrogen and carbon dioxide as cushion gas for underground gas storage reservoir. Geosyst. Eng. 2015, 18 (3), 163−167. (26) Oldenburg, C. M. Carbon Dioxide as Cushion Gas for Natural Gas Storage. Energy Fuels 2003, 17 (1), 240−246. (27) Kü hn, M.; Nakaten, N.; Streibel, M.; Kempka, T. CO2 Geological Storage and Utilization for a Carbon Neutral “Power-togas-to-power” Cycle to Even Out Fluctuations of Renewable Energy Provision. Energy Procedia 2014, 63, 8044−8049. (28) Kühn, M.; Streibel, M.; Nakaten, N.; Kempka, T. Integrated Underground Gas Storage of CO2 and CH4 to Decarbonise the “Power-to-gas-to-gas-to-power” Technology. Energy Procedia 2014, 59 (5), 9−15. (29) Niu, C.; Tan, Y. Numerical simulation and analysis of migration law of gas mixture using carbon dioxide as cushion gas in underground gas storage reservoir. J. Harbin Inst. Technol. (New Ser.) 2014, 21 (3), 121−128. (in Chinese with English abstract). (30) COMSOL Multiphysics, version 5.2a; COMSOL AB: Stockholm, 2018. www.comsol.com. (31) Wang, L.; Yu, Q. Methane adsorption on porous nano-silica in the presence of water: An experimental and ab initio study. J. Colloid Interface Sci. 2016, 467, 60−69. (32) Setzmann, U.; Wagner, W. A New Equation of State and Tables of Thermodynamic Properties for Methane Covering the Range from the Melting Line to 625 K at Pressures up to 100 MPa. J. Phys. Chem. Ref. Data 1991, 20 (6), 1061−1155. (33) Span, R.; Wagner, W. A New Equation of State for Carbon Dioxide Covering the Fluid Region from the Triple-Point Temperature to 1100 K at Pressures up to 800 MPa. J. Phys. Chem. Ref. Data 1996, 25 (6), 1509−1596. (34) Oldenburg, C.; Pruess, K.; Benson, S. M. Process modeling of CO2 injection into natural gas reservoirs for carbon sequestration and enhanced gas recovery. Energy Fuels 2001, 15 (2), 293−298. (35) Kunz, O.; Wagner, W. The GERG-2008 Wide-Range Equation of State for Natural Gases and Other Mixtures: An Expansion of GERG-2004. J. Chem. Eng. Data 2012, 57 (11), 3032−3091. (36) Friend, D. G. NIST 14. NIST Mixture Property Database. https://www.nist.gov/publications/nist-14-nist-mixture-propertydatabase. (37) Li, Q.; Ito, K.; Wu, Z.; Lowry, C. S.; Loheide, S. P., II COMSOL Multiphysics: A novel approach to ground water modeling. Groundwater 2009, 47 (4), 480−487. (38) Panfilov, M. Physicochemical Fluid Dynamics in Porous Media: Applications in Geosciences and Petroleum Engineering; John Wiley & Sons: 2018. (39) Ran, Q. Dynamic model of physical parameters in numerical reservoir simulations on hydro-mechanical coupling. Pet. Explor. Dev. 1997, 24 (3), 61−65. (40) Oldenburg, C. M.; Moridis, G. J.; Spycher, N.; Pruess, K. EOS7C Version 1.0: TOUGH2 Module for Carbon Dioxide or Nitrogen in Natural Gas (Methane) Reservoirs; Lawrence Berkeley National Laboratory, Berkeley, CA, 2004. (41) Winn, E. B. The temperature dependence of the self-diffusion coefficients of argon, neon, nitrogen, oxygen, carbon dioxide, and methane. Phys. Rev. 1950, 80 (6), 1024.

6541

DOI: 10.1021/acs.energyfuels.9b00518 Energy Fuels 2019, 33, 6527−6541