Co-Optimization of CO 2 Storage and Enhanced Gas Recovery Using Carbonated Water and Supercritical CO 2

: CO 2 -based enhanced gas recovery (EGR) is an appealing method with the dual beneﬁt of improving recovery from mature gas reservoirs and storing CO 2 in the subsurface, thereby reducing net emissions. However, CO 2 injection for EGR has the drawback of excessive mixing with the methane gas, therefore, reducing the quality of gas produced and leading to an early breakthrough of CO 2 . Although this issue has been identiﬁed as a major obstacle in CO 2 -based EGR, few strategies have been suggested to mitigate this problem. We propose a novel hybrid EGR method that involves the injection of a slug of carbonated water before beginning CO 2 injection. While still ensuring CO 2 storage, carbonated water hinders CO 2 -methane mixing and reduces CO 2 mobility, therefore delaying breakthrough. We use reservoir simulation to assess the feasibility and beneﬁt of the proposed method. Through a structured design of experiments (DoE) framework, we perform sensitivity analysis, uncertainty assessment, and optimization to identify the ideal operation and transition conditions. Results show that the proposed method only requires a small amount of carbonated water injected up to 3% pore volumes. This EGR scheme is mainly inﬂuenced by the heterogeneity of the reservoir, slug volume injected, and production rates. Through Monte Carlo simulations, we demonstrate that high recovery factors and storage ratios can be achieved while keeping recycled CO 2 ratios low.


Introduction
The global demand for energy has been rising over the past few decades and is expected to continue growing. In 2021 alone, this demand is projected to increase by approximately 4.6% [1]. Expanding industries and the world's continuously growing population are major contributors to the global thirst for energy [2]. Despite significant advancements in renewable energy technologies, fossil fuels still remain the dominant source of world energy, and this trend is predicted to continue [3][4][5]. Currently, oil, natural gas, and coal contribute more than 70% of primary energy consumed globally [6]. However, a preference for natural gas is becoming evident, as it is considered a relatively cleaner source of energy. The share of natural gas as a source of global primary energy has risen from approximately 21% to 25% over the last decade, while comparatively, the shares of oil and coal have dropped [7]. This suggests that natural gas is expected to play a significant role in the global energy mix, especially for power generation and other important processes such as water desalination [8][9][10].
In contrast to the rise in demand for natural gas, the discovery of new conventional gas reservoirs has been decreasing. The global reserves-to-production ratio of natural gas has gradually declined over the last two decades [7]. Easily recoverable gas has become more difficult to find, yet the demand for it continues to rise. Therefore, it is crucial to improve the recovery from existing reservoirs, which can be achieved through enhanced gas recovery (EGR) processes [11,12]. and displacement efficiency, and its mobility is inhibited due to the presence of carbonated water, therefore delaying breakthrough and reducing recycling.
In this work, we demonstrate the effectiveness of the proposed CO 2 -based EGR method using reservoir simulation on a synthetic case and a real gas field. Furthermore, we perform a rigorous sensitivity analysis and uncertainty assessment to identify important uncertain and operation parameters that influence the processes under study. We also carry out optimization to identify the ideal operation and transition conditions considering methane recovery factor, CO 2 stored, and CO 2 recycled as our objective functions. These objective functions are coupled in a global objective function which is maximized to cooptimize between methane recovery and CO 2 storage. The sensitivity analysis, uncertainty assessment, and optimization are carried out as part of a structured Design of Experiments (DoE) framework that approaches the problem in a non-deterministic manner.
This paper is organized as follows. We first review the governing equations of multiphase fluid flow in porous media. The reservoir models are then introduced, followed by the proposed modeling workflow using a DoE. In the Results section, the thermodynamics of CO 2 /water systems corresponding to CO 2 solubility and phase behavior is validated with measured data. The optimization approach is then applied for a synthetic case and a 3D field case. Finally, we provide the main conclusions.

Flow Model
The governing equations for multiphase compositional flow in porous media are given by the component balance equations, Darcy's law, and the thermodynamic equilibrium between the phases [36]. In this work, we neglect capillary pressure due to the negligible interfacial tension in CO 2 -methane systems at the conditions of interest in this study [37].
The material balance for each component in the phases is given by: where φ is the porosity, c is the overall molar density, z i is the overall mole fraction of component i, U i is the molar flux of component i, and n c is the total number of components. The molar flux U i is given by: where α is the phase (water, oil, or gas), c α is the molar density of phase α, x i,α is the mole fraction of component i in phase α, v α is the phase velocity, S α is the saturation of phase α, and J i,α is the diffusion flux of component i in phase α. The global and phase compositions z i and x i,α respectively are constrained by the relation: The velocity of each phase is given by Darcy's law [38] as: The equation for pressure is written based on the concept of volume balance [41]: where ξ is the total fluid compressibility and v i is the total partial molar volume of component i. The local thermodynamic equilibrium is defined by the equality of the fugacities of each component in the phases α 1 and α 2 , that is: f i,α 1 (T, p, x j,α 1 ) = f i,α 2 (T, p, x j,α 2 ), i = 1, . . . , n c , j = 1, . . . , n c − 1 (9) In the above equation, f i,α is the fugacity of a component i in phase α, which can be calculated using an equation of state.

CO 2 -Water Phase Behavior
Understanding the interaction of CO 2 with water and the thermodynamic properties of the resulting mixture is an important part of CO 2 storage research. The dissolution of CO 2 in formation water is a primary CO 2 trapping mechanism over the long term. This reaction is relatively fast, and therefore, thermodynamic equilibrium is assumed. This equilibrium is expressed as the equality of fugacities in the gas and aqueous phases, such that: where, f i,g is the fugacity of the ith component in the gas phase, and f i,aq is the fugacity of the ith component in the aqueous phase. For CO 2 gas, its fugacity in the gas phase is related to its partial pressure by: where ϕ is the fugacity coefficient calculated from the Peng-Robinson EOS [42], and P CO 2 is the partial pressure of CO 2 . In the aqueous phase, the fugacity of CO 2 is calculated using extended Henry's law [43].
where, y CO 2 is the mole fraction of CO 2 in the aqueous phase and H CO 2 is Henry's constant. In pure water, Henry's constant in Equation (12) is calculated as a function of pressure, temperature, the universal gas constant, and partial molar volume with respect to a reference pressure, is given by: In saline water, a correction is made to Henry's constant to account for salinity, as follows: where H s,CO 2 is Henry's constant for CO 2 in the brine, H CO 2 is Henry's constant for CO 2 in pure water, k s,CO 2 is the Setchenov salting-out coefficient for a CO 2 -water aqueous phase, and m s is the molality of the dissolved salt.

Mechanistic Simulations
In this study, we use CMG GEM v2019.1, a commercial simulator from Computer Modeling Group (CMG). GEM is a fully-coupled compositional Equation of State (EOS) simulator capable of modeling subsurface flow problems, including CO 2 storage in oil and gas reservoirs [44].

Synthetic Reservoir 2D Cross-Section
We first consider a 2D gas reservoir model of physical dimensions 500 × 10 × 100 m representing a vertical (I-K) cross-section. The cross-section is discretized into a 50 × 1 × 20 regular Cartesian grid, as shown in Figure 1. There are 10 geological layers with varying permeability. Each layer is assumed to be homogenous, and the physical properties are uniform. The heterogeneity of the reservoir model is quantified by the Dykstra-Parsons coefficient [45], which is 0.5 for the case shown. There are two vertical wells (an injector and producer) completed across the entire thickness of the reservoir. Other physical, initialization, and sensitivity parameters for the simulation model are summarized in Appendix A, Table A1. The simulations of CO 2 -EGR and storage for the synthetic reservoir cross-section are controlled by the constraints summarized in Table A2.
In saline water, a correction is made to Henry's constant to account for salinity, as follows:

Mechanistic Simulations
In this study, we use CMG GEM v2019.1, a commercial simulator from Computer Modeling Group (CMG). GEM is a fully-coupled compositional Equation of State (EOS) simulator capable of modeling subsurface flow problems, including CO2 storage in oil and gas reservoirs [44].

Synthetic Reservoir 2D Cross-Section
We first consider a 2D gas reservoir model of physical dimensions 500 10 100 m × × representing a vertical (I-K) cross-section. The cross-section is discretized into a 50 1 20 × × regular Cartesian grid, as shown in Figure 1. There are 10 geological layers with varying permeability. Each layer is assumed to be homogenous, and the physical properties are uniform. The heterogeneity of the reservoir model is quantified by the Dykstra-Parsons coefficient [45], which is 0.5 for the case shown. There are two vertical wells (an injector and producer) completed across the entire thickness of the reservoir. Other physical, initialization, and sensitivity parameters for the simulation model are summarized in Appendix A, Table A1. The simulations of CO2-EGR and storage for the synthetic reservoir cross-section are controlled by the constraints summarized in Table  A2. For consistency, we define two dimensionless ratios that capture the production and injection rates and the pressure at which EGR has commenced. These ratios are defined as: • Productivity ratio-is the ratio of the production rate to the injection rate. • Depletion pressure ratio-is the ratio of the pressure at the start of EGR to the original reservoir pressure.

Real Reservoir 3D Sector
To demonstrate the proposed method in the presence of complex heterogeneities, we consider a 3D sector from a real gas field as shown in Figure 2. The sector is approximately 3500 × 4100 × 500 ft in dimension and is discretized into a 63 × 75 × 32 Cartesian grid, resulting in 151,200 grid blocks. The porosity varies between 0.05 and 0.23, while the permeability varies from less than 1 mD to more than 1600 mD. There is one vertical producer and one vertical injector. The producer is placed up-dip from the location of the injector. Relevant physical and initialization parameters are summarized in Table A3. Simulation constraints are presented in Table A4. consider a 3D sector from a real gas field as shown in Figure 2. The sector is approximately 3500 4100 500 ft × × in dimension and is discretized into a 63 75 32 × × Cartesian grid, resulting in 151,200 grid blocks. The porosity varies between 0.05 and 0.23, while the permeability varies from less than 1 mD to more than 1600 mD. There is one vertical producer and one vertical injector. The producer is placed up-dip from the location of the injector. Relevant physical and initialization parameters are summarized in Table A3. Simulation constraints are presented in Table A4.

DoE Framework
We defined a Design of Experiments (DoE) workflow to perform sensitivity analysis, uncertainty assessment, and optimization. The key steps of this workflow are presented in Figure 3. We began by identifying and screening multiple possible uncertainty parameters. We used 2-level Plackett-Burman analysis [46,47] because it is computationally inexpensive [48]. However, in some cases, this design might not be enough to capture the response surface accurately. Therefore, we switched to Latin Hypercube Sampling (LHS). Using the experiments generated through LHS, we performed Sobol analysis [49] to identify the most sensitive parameters. These are the parameters considered for the uncertainty assessment. Furthermore, using the results of LHS sampling, we built a polynomial proxy model. We verified the quality of the proxy model using blind tests. Once the quality of the proxy model was verified, Monte Carlo simulations were performed using the proxy as part of the uncertainty assessment. We also performed optimization using the proxy model as part of the Latin Hypercube plus Proxy Optimization algorithm. The flowchart of the optimization algorithm is presented in Figure 4.

DoE Framework
We defined a Design of Experiments (DoE) workflow to perform sensitivity analysis, uncertainty assessment, and optimization. The key steps of this workflow are presented in Figure 3. We began by identifying and screening multiple possible uncertainty parameters. We used 2-level Plackett-Burman analysis [46,47] because it is computationally inexpensive [48]. However, in some cases, this design might not be enough to capture the response surface accurately. Therefore, we switched to Latin Hypercube Sampling (LHS). Using the experiments generated through LHS, we performed Sobol analysis [49] to identify the most sensitive parameters. These are the parameters considered for the uncertainty assessment. Furthermore, using the results of LHS sampling, we built a polynomial proxy model. We verified the quality of the proxy model using blind tests. Once the quality of the proxy model was verified, Monte Carlo simulations were performed using the proxy as part of the uncertainty assessment. We also performed optimization using the proxy model as part of the Latin Hypercube plus Proxy Optimization algorithm. The flowchart of the optimization algorithm is presented in Figure 4.
The proposed EGR method was evaluated based on 3 main dimensionless objective functions; methane recovery factor, CO 2 stored, and CO 2 recycled. These objective functions are defined as: Recovery Factor = CH 4 mass produced CH 4 mass originally in place (15) CO 2 Stored = CO 2 mass injected − CO 2 mass produced Maximum storable CO 2 mass at target p, T CO 2 Recycled = CO 2 mass produced CO 2 mass injected (17) These objective functions are coupled in a linear global objective function of the form aX 1 + bX 2 − cX 3 , which is maximized to identify the ideal operation and transition conditions. The variables,X 1 , X 2 , and X 3 are the recovery factor, CO 2 stored, and CO 2 recycled respectively while a, b, and c are bios coefficients that could be varied to reflect the relative importance of the objective functions in the optimization. aX bX cX + − , which is maximized to identify the ideal operation and transition conditions. The variables, 1 X , 2 X , and 3 X are the recovery factor, CO2 stored, and CO2 recycled respectively while a , b , and c are bios coefficients that could be varied to reflect the relative importance of the objective functions in the optimization.

Simulator Validation
Prior to conducting multi-dimensional simulations, we performed controlled calibration and verification of the simulator used to identify and assess the key governing mechanism related to CO2-EGR and storage. These mechanisms were studied to ensure the accuracy and representativeness of our simulation models. The investigated mechanisms include: -CO2 solubility in water as a function of pressure, temperature, and salinity. -CO2 density (molar volume) as a function of pressure and temperature. -CO2-saturated water density as a function of pressure and temperature. -Water vaporization as a function of pressure and temperature.
The solubility of CO2 in pure water versus pressure and temperature is shown in Figure 5. CO2 solubility increases with increasing pressure and decreases with increasing temperature. The simulator calculates the fugacity of gas components soluble in the aqueous phase using Henry's law [43], as previously discussed. There is a good match between the experimental data and the calculations for the range of data presented. The model is also capable of calculating CO2 solubility in water at different salinities (see Figure 6),

Simulator Validation
Prior to conducting multi-dimensional simulations, we performed controlled calibration and verification of the simulator used to identify and assess the key governing mechanism related to CO 2 -EGR and storage. These mechanisms were studied to ensure the accuracy and representativeness of our simulation models. The investigated mechanisms include: -CO 2 solubility in water as a function of pressure, temperature, and salinity. The solubility of CO 2 in pure water versus pressure and temperature is shown in Figure 5. CO 2 solubility increases with increasing pressure and decreases with increasing temperature. The simulator calculates the fugacity of gas components soluble in the aqueous phase using Henry's law [43], as previously discussed. There is a good match between the experimental data and the calculations for the range of data presented. The model is also capable of calculating CO 2 solubility in water at different salinities (see Figure 6), where the method based on the Scaled Particle Theory by Li and Nghiem [51] is used to correct Henry's constants in the presence of salt.
Energies 2021, 14, x FOR PEER REVIEW 9 of 23 Figure 5. Solubility of CO2 in pure water as a function of pressure and temperature. The discrete points represent experimental data from [52], and the continuous lines are calculated by the simulator. Figure 6. Solubility of CO2 in 0, 1, and 3 molal brine at 50 °C (left) and 100 °C (right) as a function of pressure. The discrete points are experimental data presented in [53], and the continuous lines are calculated by the simulator.
The density of the gaseous phase is calculated using the Peng-Robinson correlation [42]. Figure 7 shows excellent agreement between the calculated density of CO2, related to the molar volume by molar mass, and the experimental data for a range of pressures and temperatures. The density of the aqueous phase considering dissolved components is calculated using the Rowe-Chou correlation [54]. Figure 8 shows the density of the CO2-saturated aqueous phase as a function of temperature and pressure. The aqueous density increases with increasing pressure and decreases with increasing temperature.
Water vaporization enables the mobilization of previously immobile water at low saturations, which could lead to salt precipitation [55,56]. The water content in the CO2rich gas phase as a function of temperature and pressure is shown in Figure 9, with a good match between the calculations and the experimental data. Figure 5. Solubility of CO 2 in pure water as a function of pressure and temperature. The discrete points represent experimental data from [52], and the continuous lines are calculated by the simulator. . Solubility of CO2 in pure water as a function of pressure and temperature. The discrete points represent experimental data from [52], and the continuous lines are calculated by the simulator. Figure 6. Solubility of CO2 in 0, 1, and 3 molal brine at 50 °C (left) and 100 °C (right) as a function of pressure. The discrete points are experimental data presented in [53], and the continuous lines are calculated by the simulator.
The density of the gaseous phase is calculated using the Peng-Robinson correlation [42]. Figure 7 shows excellent agreement between the calculated density of CO2, related to the molar volume by molar mass, and the experimental data for a range of pressures and temperatures. The density of the aqueous phase considering dissolved components is calculated using the Rowe-Chou correlation [54]. Figure 8 shows the density of the CO2-saturated aqueous phase as a function of temperature and pressure. The aqueous density increases with increasing pressure and decreases with increasing temperature.
Water vaporization enables the mobilization of previously immobile water at low saturations, which could lead to salt precipitation [55,56]. The water content in the CO2rich gas phase as a function of temperature and pressure is shown in Figure 9, with a good match between the calculations and the experimental data. Figure 6. Solubility of CO 2 in 0, 1, and 3 molal brine at 50 • C (left) and 100 • C (right) as a function of pressure. The discrete points are experimental data presented in [53], and the continuous lines are calculated by the simulator.
The density of the gaseous phase is calculated using the Peng-Robinson correlation [42]. Figure 7 shows excellent agreement between the calculated density of CO 2 , related to the molar volume by molar mass, and the experimental data for a range of pressures and temperatures. The density of the aqueous phase considering dissolved components is calculated using the Rowe-Chou correlation [54]. Figure 8 shows the density of the CO 2 -saturated aqueous phase as a function of temperature and pressure. The aqueous density increases with increasing pressure and decreases with increasing temperature.
Water vaporization enables the mobilization of previously immobile water at low saturations, which could lead to salt precipitation [55,56]. The water content in the CO 2 -rich gas phase as a function of temperature and pressure is shown in Figure 9, with a good match between the calculations and the experimental data. temperatures. The density of the aqueous phase considering dissolved components is cal-culated using the Rowe-Chou correlation [54]. Figure 8 shows the density of the CO2-saturated aqueous phase as a function of temperature and pressure. The aqueous density increases with increasing pressure and decreases with increasing temperature.
Water vaporization enables the mobilization of previously immobile water at low saturations, which could lead to salt precipitation [55,56]. The water content in the CO2rich gas phase as a function of temperature and pressure is shown in Figure 9, with a good match between the calculations and the experimental data.  . Molar volume of CO2 as a function of temperature and pressure. The discrete points are experimental data presented in [52], and the continuous lines are calculated by the simulator.

Figure 8.
Density of CO2-saturated water as a function of temperature and pressure. The discrete points are experimental data presented in [57], and the continuous lines are calculated by the simulator.

Figure 9.
Mole fraction of water in the CO2-rich gas phase in a water-CO2 mixture as a function of temperature and pressure. The discrete points are experimental data from [52], and the continuous lines are calculated by the simulator.

Benefit of Carbonated Water Injection
We performed mechanistic simulations to assess the feasibility of the proposed EGR method, including sensitivity analysis to quantify the uncertainty and identify the optimum transition and operation conditions. Each simulation covers four main periods. First, natural depletion followed by carbonated water injection when the reservoir pressure drops to a certain depletion pressure limit. After a defined volume of carbonated water is injected, pure CO2 injection follows. Once the reservoir pressure reaches the original reservoir pressure, injection is stopped.
Before carrying out a detailed sensitivity analysis, uncertainty assessment, and optimization, we compared two baseline cases to observe the difference between EGR that begins with pure CO2 injection and injection of a slug of carbonated water first. We observe the flow patterns in the 2D cross-sections comparing the two cases, as shown in Figure 10.   [57], and the continuous lines are calculated by the simulator. Figure 9. Mole fraction of water in the CO2-rich gas phase in a water-CO2 mixture as a function of temperature and pressure. The discrete points are experimental data from [52], and the continuous lines are calculated by the simulator.

Benefit of Carbonated Water Injection
We performed mechanistic simulations to assess the feasibility of the proposed EGR method, including sensitivity analysis to quantify the uncertainty and identify the optimum transition and operation conditions. Each simulation covers four main periods. First, natural depletion followed by carbonated water injection when the reservoir pressure drops to a certain depletion pressure limit. After a defined volume of carbonated water is injected, pure CO2 injection follows. Once the reservoir pressure reaches the original reservoir pressure, injection is stopped.
Before carrying out a detailed sensitivity analysis, uncertainty assessment, and optimization, we compared two baseline cases to observe the difference between EGR that begins with pure CO2 injection and injection of a slug of carbonated water first. We observe the flow patterns in the 2D cross-sections comparing the two cases, as shown in Figure 10. . Mole fraction of water in the CO 2 -rich gas phase in a water-CO 2 mixture as a function of temperature and pressure. The discrete points are experimental data from [52], and the continuous lines are calculated by the simulator.

Benefit of Carbonated Water Injection
We performed mechanistic simulations to assess the feasibility of the proposed EGR method, including sensitivity analysis to quantify the uncertainty and identify the optimum transition and operation conditions. Each simulation covers four main periods. First, natural depletion followed by carbonated water injection when the reservoir pressure drops to a certain depletion pressure limit. After a defined volume of carbonated water is injected, pure CO 2 injection follows. Once the reservoir pressure reaches the original reservoir pressure, injection is stopped.
Before carrying out a detailed sensitivity analysis, uncertainty assessment, and optimization, we compared two baseline cases to observe the difference between EGR that begins with pure CO 2 injection and injection of a slug of carbonated water first. We observe the flow patterns in the 2D cross-sections comparing the two cases, as shown in Figure 10.  Figure 10 shows a series of concentration maps comparing the flow and sweep patterns of the case of CO2 injection only and the hybrid case of carbonated water combined with CO2 injection. CO2 injection and carbonated water injection begin in their respective cases. After 5 months, better sweep and displacement efficiency can be seen in the hybrid case. In the hybrid case, the high permeability channels contain water which inhibits the mobility of CO2 through these channels, unlike in the absence of carbonated water injection. After 13 months, the injection of the pre-determined carbonated water slug is completed, and the hybrid case switches to pure CO2 injection. At 61 months, CO2 breaks through at the production well in the case with no carbonated water slug injected. Comparatively, CO2 breaks through after 77 months in the hybrid case.
To understand the effect of the volume of carbonated water injected, we performed several simulations with different slug volumes injected from 0 up to 0.1 pore volumes.  Figure 10 shows a series of concentration maps comparing the flow and sweep patterns of the case of CO 2 injection only and the hybrid case of carbonated water combined with CO 2 injection. CO 2 injection and carbonated water injection begin in their respective cases. After 5 months, better sweep and displacement efficiency can be seen in the hybrid case. In the hybrid case, the high permeability channels contain water which inhibits the mobility of CO 2 through these channels, unlike in the absence of carbonated water injection. After 13 months, the injection of the pre-determined carbonated water slug is completed, and the hybrid case switches to pure CO 2 injection. At 61 months, CO 2 breaks through at the production well in the case with no carbonated water slug injected. Comparatively, CO 2 breaks through after 77 months in the hybrid case.
To understand the effect of the volume of carbonated water injected, we performed several simulations with different slug volumes injected from 0 up to 0.1 pore volumes. We compare these cases by the recovery factor, CO 2 stored, and the number of pore volumes injected at the breakthrough time, which is a direct dimensionless indicator of breakthrough time, as shown in Figure 11. The recovery factor slightly improves as more carbonated water is injected. This behavior is attributed to the lower mobility of water, thus giving a better displacement efficiency. However, the CO 2 stored decreases if more carbonated water is injected, which is attributed to carbonated water occupying potential storage space in the reservoir. On the other hand, breakthrough is delayed for up to 0.1 PV of carbonated water injected. The delay, however, approaches a plateau as more carbonated water is injected. This points to a decreasing benefit as more water is injected. We compare these cases by the recovery factor, CO2 stored, and the number of pore volumes injected at the breakthrough time, which is a direct dimensionless indicator of breakthrough time, as shown in Figure 11. The recovery factor slightly improves as more carbonated water is injected. This behavior is attributed to the lower mobility of water, thus giving a better displacement efficiency. However, the CO2 stored decreases if more carbonated water is injected, which is attributed to carbonated water occupying potential storage space in the reservoir. On the other hand, breakthrough is delayed for up to 0.1 PV of carbonated water injected. The delay, however, approaches a plateau as more carbonated water is injected. This points to a decreasing benefit as more water is injected. Considering the opposing trends observed in Figure 11, we couple these phenomena in a single objective function to quantify the overall benefit of the proposed hybrid scheme. When this objective function is plotted as a function of carbonated water slug volume, we observe an increasing trend followed by a decrease, as shown in Figure 12. This suggests that for up to some amount of carbonated water injected (~0.3 PV in this case), there is a benefit. However, after this amount, there is no evident benefit in injecting a slug of carbonated water considering recovery factor, CO2 stored, and breakthrough time. Considering the opposing trends observed in Figure 11, we couple these phenomena in a single objective function to quantify the overall benefit of the proposed hybrid scheme. When this objective function is plotted as a function of carbonated water slug volume, we observe an increasing trend followed by a decrease, as shown in Figure 12. This suggests that for up to some amount of carbonated water injected (~0.3 PV in this case), there is a benefit. However, after this amount, there is no evident benefit in injecting a slug of carbonated water considering recovery factor, CO 2 stored, and breakthrough time.

Sensitivity Analysis, Uncertainty Assessment, and Optimization
To understand the most dominant parameters influencing the proposed EGR method, we initially performed sensitivity analysis, uncertainty assessment, and optimization on the synthetic cross-section. This study was done following the DoE framework presented in Section 4. The parameters analyzed are presented in Table A1. As part of the

Sensitivity Analysis, Uncertainty Assessment, and Optimization
To understand the most dominant parameters influencing the proposed EGR method, we initially performed sensitivity analysis, uncertainty assessment, and optimization on the synthetic cross-section. This study was done following the DoE framework presented in Section 4. The parameters analyzed are presented in Table A1. As part of the sensitivity analysis, the influence of 9 parameters on the methane recovery factor, CO 2 stored, and CO 2 recycled was investigated. The sensitivity was quantified by the Sobol method [49,58], which is a variance-based sensitivity analysis method. About 108 experiments sampled using Latin Hypercube sampling were performed and were sufficient to capture the response surface. The results of the sensitivity analysis are presented in Figure 13, where DP is the Dykstra-Parsons coefficient to quantify heterogeneity, and "EGR Comm. Pressure" is the pressure at which EGR begins.

Sensitivity Analysis, Uncertainty Assessment, and Optimization
To understand the most dominant parameters influencing the proposed EGR method, we initially performed sensitivity analysis, uncertainty assessment, and optimization on the synthetic cross-section. This study was done following the DoE framework presented in Section 4. The parameters analyzed are presented in Table A1. As part of the sensitivity analysis, the influence of 9 parameters on the methane recovery factor, CO2 stored, and CO2 recycled was investigated. The sensitivity was quantified by the Sobol method [49,58], which is a variance-based sensitivity analysis method. About 108 experiments sampled using Latin Hypercube sampling were performed and were sufficient to capture the response surface. The results of the sensitivity analysis are presented in Figure  13, where DP is the Dykstra-Parsons coefficient to quantify heterogeneity, and "EGR Comm. Pressure" is the pressure at which EGR begins.  A 1% significance threshold was set to eliminate parameters, which did not influence the objective functions. Considering this criterion, parameters chosen for uncertainty assessment are: -Dykstra-Parsons coefficient; -Production rate; -EGR commencement pressure; -Carbonated water slug volume injected; -Average permeability; -Vertical to horizontal permeability ratio.
In the uncertainty quantification, a polynomial proxy model was built using 300 training experiments and verified using 14 blind test cases. Figure 14 shows the proxy verification for the main objective functions; recovery factor, CO 2 stored, and CO 2 recycled. The proxies generated were accurate enough to be used in Monte Carlo simulations, with R 2 values above 0.9.
In the uncertainty quantification, a polynomial proxy model was built using 300 training experiments and verified using 14 blind test cases. Figure 14 shows the proxy verification for the main objective functions; recovery factor, CO2 stored, and CO2 recycled. The proxies generated were accurate enough to be used in Monte Carlo simulations, with R 2 values above 0.9. In the next step, the verified proxy was used for Monte Carlo (MC) simulation. Each MC simulation uses 65,000 samples to ensure few gaps in the sampling space. Better space filling in the sampling space generally gives better predictions [59,60]. The results for methane recovery factor, CO2 stored, and CO2 recycled are shown in Figure 15. From the MC simulation, the most likely methane recovery factor is ~94%. This is in agreement with recovery factors from some gas reservoirs around the world, which have reached up to ~96% [61]. The possible recovery factors range between 90% and 98%. The most likely value for CO2 stored as a percent of the maximum storable CO2 at the reservoir conditions is ~79%. The range of possible storage ratios is from 71% to 87%. On the other hand, the most likely value for CO2 recycled is ~16% of the injected CO2, while the possible range is from 10% to 30%. These results are promising as they show a likelihood of achieving high recovery factors while storing most of the CO2 injected using the proposed EGR method. In the next step, the verified proxy was used for Monte Carlo (MC) simulation. Each MC simulation uses 65,000 samples to ensure few gaps in the sampling space. Better space filling in the sampling space generally gives better predictions [59,60]. The results for methane recovery factor, CO 2 stored, and CO 2 recycled are shown in Figure 15. From the MC simulation, the most likely methane recovery factor is~94%. This is in agreement with recovery factors from some gas reservoirs around the world, which have reached up tõ 96% [61]. The possible recovery factors range between 90% and 98%. The most likely value for CO 2 stored as a percent of the maximum storable CO 2 at the reservoir conditions is~79%. The range of possible storage ratios is from 71% to 87%. On the other hand, the most likely value for CO 2 recycled is~16% of the injected CO 2 , while the possible range is from 10% to 30%. These results are promising as they show a likelihood of achieving high recovery factors while storing most of the CO 2 injected using the proposed EGR method. The polynomial proxy model built was used in the optimization algorithm described in Section 4 to identify the optimum operation and transition conditions for the proposed EGR method. This process is done by maximizing a global objective function that couples the recovery factor, CO2 stored, and CO2 recycled. The algorithm analyzes 500 experiments within the optimization framework, as shown in Figure 16.  The polynomial proxy model built was used in the optimization algorithm described in Section 4 to identify the optimum operation and transition conditions for the proposed EGR method. This process is done by maximizing a global objective function that couples the recovery factor, CO 2 stored, and CO 2 recycled. The algorithm analyzes 500 experiments within the optimization framework, as shown in Figure 16. The polynomial proxy model built was used in the optimization algorithm described in Section 4 to identify the optimum operation and transition conditions for the proposed EGR method. This process is done by maximizing a global objective function that couples the recovery factor, CO2 stored, and CO2 recycled. The algorithm analyzes 500 experiments within the optimization framework, as shown in Figure 16. Initially, a set of experiments is generated by Latin Hypercube sampling and is used to train a proxy model. Once the proxy model is generated, it tries to find an optimum solution, which is validated by a full-physics simulation. In Figure 16, there is a jump in solution quality when the proxy model begins the optimization at ~300 experiments. This result shows the efficiency of the proxy model in finding possible optimum solutions. Some outlier points are the result of the optimizer searching possible far-field solutions to ensure it is not stuck at local maxima. Initially, a set of experiments is generated by Latin Hypercube sampling and is used to train a proxy model. Once the proxy model is generated, it tries to find an optimum solution, which is validated by a full-physics simulation. In Figure 16, there is a jump in solution quality when the proxy model begins the optimization at~300 experiments. This result shows the efficiency of the proxy model in finding possible optimum solutions. Some outlier points are the result of the optimizer searching possible far-field solutions to ensure it is not stuck at local maxima.
From the optimization study, the operation and transition conditions of the optimum case were: -Productivity ratio: 0.5; -Depletion pressure ratio: 0.115; -Carbonated water PV injected: 1.5%.
Comparing cases close to the optimum case, we identified that the lower productivity and depletion pressure ratios resulted in the highest global objective function values. In addition, a small PV ranging from 1% to 3% of carbonated water was beneficial. These conditions resulted in the delayed breakthrough of CO 2 and relatively higher methane recovery factors.

Field Case (Viking Field)
Viking Field is a gas field off the coast of Lincolnshire, England, in the North Sea [62]. It lies in the Rotliegendes sandstone at an approximate depth of 9000 ft subsea [63]. Viking Field has two main reservoirs which are currently depleted, Viking A and Viking B, and several smaller gas pools which totally held approximately 3.6 Tcf of initial gas in place. The produced gas from Viking Field contains 0.3-0.35 mol% of CO 2 [64].
To demonstrate the proposed hybrid EGR method on a more complex multi-dimensional domain, two cases were compared on a sector from the Viking B field. For computational efficiency, the depletion period was omitted, and the initialization reservoir pressure was set to 300 psi, reflecting the actual abundant pressure of the depleted gas reservoir. Figure 17 shows the flow behavior in the sector as CO 2 displaces methane and moves from the injection well to the production well. The wells are placed in a quarter 5-spot pattern. At time zero, injection begins in the respective cases. At 6 months, it is clear how the more mobile CO 2 is advancing faster through the higher permeability layers compared to the carbonated water. At 22 months, carbonated water injection is completed, and pure CO 2 injection begins in the hybrid case. At 129 months, breakthrough occurs in the CO 2 -only case, and the production well is shut-in due to excessive CO 2 production. In comparison, CO 2 breaks through in the hybrid case after 208 months, and the production well is also shut-in.
sure was set to 300 psi, reflecting the actual abundant pressure of the depleted gas reservoir. Figure 17 shows the flow behavior in the sector as CO2 displaces methane and moves from the injection well to the production well. The wells are placed in a quarter 5-spot pattern. At time zero, injection begins in the respective cases. At 6 months, it is clear how the more mobile CO2 is advancing faster through the higher permeability layers compared to the carbonated water. At 22 months, carbonated water injection is completed, and pure CO2 injection begins in the hybrid case. At 129 months, breakthrough occurs in the CO2only case, and the production well is shut-in due to excessive CO2 production. In comparison, CO2 breaks through in the hybrid case after 208 months, and the production well is also shut-in.  The above demonstration of flow behavior and breakthrough time on a sector from the Viking Field shows the potential applicability of the proposed hybrid method on more complex multi-dimensional domains. It is, however, important to identify some of the uncertainties associated with scaling up to the field level. The heterogeneity at the field scale is a major source of uncertainty and directly impacts the flow behavior and breakthrough of injected fluids. The available pore volume is also an inherent source of uncertainty, as it affects the available space for CO2 storage. Furthermore, certain operational and decisional parameters, for instance, well spacing, affect breakthrough and are parameters to be considered in order to optimize the implementation of the method.

Potential Field Applications, Economics, and Limitations
The proposed hybrid method is very novel and has not been deployed in any actual gas fields as of now. Furthermore, this method has not been previously suggested in published literature, to the best of our knowledge. However, co-injection of carbonated water with supercritical CO2 for EGR is quite promising, as shown by the presented simulation studies and the field demonstration using a sector from the Viking gas field. Based on Figure 17. Global CO 2 mole fraction at different times for CO 2 injection only (left) and carbonated water combined with CO 2 injection (right). The time is given in months since injection began.
The above demonstration of flow behavior and breakthrough time on a sector from the Viking Field shows the potential applicability of the proposed hybrid method on more complex multi-dimensional domains. It is, however, important to identify some of the uncertainties associated with scaling up to the field level. The heterogeneity at the field scale is a major source of uncertainty and directly impacts the flow behavior and breakthrough of injected fluids. The available pore volume is also an inherent source of uncertainty, as it affects the available space for CO 2 storage. Furthermore, certain operational and decisional parameters, for instance, well spacing, affect breakthrough and are parameters to be considered in order to optimize the implementation of the method.

Potential Field Applications, Economics, and Limitations
The proposed hybrid method is very novel and has not been deployed in any actual gas fields as of now. Furthermore, this method has not been previously suggested in published literature, to the best of our knowledge. However, co-injection of carbonated water with supercritical CO 2 for EGR is quite promising, as shown by the presented simulation studies and the field demonstration using a sector from the Viking gas field. Based on published works from the literature, we identify other candidate reservoirs which could potentially benefit from the proposed hybrid method: • Altmark Field in Germany is a gas field that has reached the tail end of its production phase, with a recovery factor > 78%. The reservoir pressure in the more productive Altensazwedel compartment of this field has dropped to 4000 kPa. Altmark Field has been considered as a candidate for CO 2 EGR coupled with storage, and several studies have been done to investigate the feasibility of this [18,65,66]. The permeability distribution in the Altensazwedel compartment is such that higher permeability layers are lower in the compartment. Studies showed that the CO 2 distribution in the reservoir after injection was strongly linked to the heterogeneity configuration, and breakthrough was likely to occur after~96 months from the beginning of injection [67]. • K12-B Field is an offshore gas field located in the Dutch sector of the North Sea, which was close to depletion when CO 2 injection field tests were conducted [68,69]. Gas production was from the Rotliegendes formation, with a heterogeneity configuration similar to the Viking Field. Monitoring studies for the CO 2 injection pilot showed CO 2 breakthrough to the production wells after 130 days for the first well and after 463 days for the second well [70]. • Otway CO 2 storage project was the first CCS project implemented in Australia. Natural gas was produced from the Naylor Field, and CO 2 injection was done into a 25-30 m thick sandstone reservoir in this field [71,72]. Monitoring studies showed the injected CO 2 form a plume that migrated up-dip to the production wells through the more permeable layers in the mid-bottom of the reservoir. Tracers injected along with the CO 2 were observed at the production wells after~150 days, pointing to an early breakthrough [73].
From a technical perspective, the obtained results support the benefit of the hybrid CO 2 injection method. However, for implementation in actual fields, the effectiveness of this method will be determined not only by its technical efficacy but rather by considering the additional benefit it provides in terms of project economics, with or without external incentives. The additional benefit given by the delayed breakthrough and extended production period must be compared against the capital and operation cost of injecting carbonated water. Various potential fields have different infrastructures in place, and the additional cost value may vary significantly depending on the already available facilities and the individual need for upgrades. Therefore, each candidate field must be evaluated individually to weigh the technical benefit of injecting carbonated water against the financial cost of its implementation.
While the proposed hybrid method in this work is promising and shows good potential, it is necessary to highlight some of the limitations or possible challenges that could be encountered while implementing this method. Some of these challenges include:

•
The proposed method is dependent on the individual reservoir heterogeneity. Therefore, each case must be individually evaluated in detail to determine whether this method could be applicable and beneficial or not.

•
Results have shown that up to 3% PV of injected carbonated water is beneficial. This is a significant amount of water required for injection and could be a challenge for fields that do not have an existing water supply available.

•
The preparation of carbonated water in large quantities under desired pressure and temperature conditions may require specialized processing facilities. This could be an additional cost to the overall project expenses.

•
The dissolution of CO 2 in water leads to the formation of carbonic acid. Carbonic acid could lead to corrosion of steel tubes that are not corrosion-resistant.

Conclusions
In this work, we propose a novel CO 2 -based enhanced gas recovery (EGR) method to mitigate the issue of early CO 2 breakthrough and excessive recycling in CO 2 EGR. We use reservoir simulation to assess the feasibility and benefit of this method. We also define a Design of Experiments (DoE) framework to understand the most dominant reservoir and operational parameters affecting the proposed method. This DoE framework is also used to quantify the uncertainty involved and perform optimization to find ideal operation and transition conditions for the proposed method. The most salient conclusions from this work include:

•
Injection of a slug of carbonated water followed by CO 2 injection is beneficial to some extent considering three factors; recovery factor, CO 2 stored, and CO 2 breakthrough/recycling. In this study, injection of up to~3% PV of carbonated water was most beneficial considering all three factors. Injection volumes more than this resulted in no benefit compared to only CO 2 injection.

•
The recovery factor was mostly affected by the reservoir heterogeneity, depletion pressure, and production rates. The CO 2 stored was mainly affected by the pore volume of carbonated water injected. This is because injection of any amount of water occupied potential storage space in the reservoir. The CO 2 breakthrough and recycling were primarily affected by the production rates and reservoir heterogeneity as high permeability paths allowed much faster CO 2 breakthrough.

•
In general, the proposed hybrid EGR method is significantly influenced by reservoir heterogeneity. Compared to only CO 2 injection, the method performs best in reservoirs with high permeability paths in the bottom layers. For CO 2 floods in gas reservoirs, CO 2 is denser than methane and flows to the bottom of the reservoir due to gravity segregation. When it encounters high permeability paths at the lower layers, it is likely to break through much faster. Injection of carbonated water becomes beneficial in such reservoirs because water, also denser than methane and the in-situ brine, flows down to the bottom of the reservoir. This water plugs the high permeability channels and thus hinders the mobility of CO 2 , therefore, delaying breakthrough. This also forces the CO 2 injected after the water to go through lower permeability layers resulting in more even sweep profiles. • From the Monte Carlo simulations performed, the recovery factor is generally high, and P10, P50, and P90 values fall above 90%. The CO 2 stored ranges between~70% to~90% of the maximum storable CO 2 at the reservoir conditions. The CO 2 recycled ranges between~10% to~30%. These results are promising and show that with the proposed EGR method, high recovery factors can be achieved while storing most of the CO 2 injected with low recycling rates. • While the proposed method is promising in terms of improving recovery with high storage ratios, it is necessary to weigh this benefit against the cost of implementing the method. This should be done considering both environmental value and external incentives in order to get an overall quantification of its effectiveness and feasibility. • Modern reservoir simulators are powerful tools that can contribute to CO 2 storage research advancement. These simulators incorporate many mathematical and empirical models that capture different physical phenomena. However, these models are limited to particular ranges of conditions. It is crucial to first verify and properly calibrate simulation models to ensure representativeness and reliability of results before moving to complex multi-dimensional problems. In this study, we take the time to properly validate the simulator used and carefully select appropriate empirical models to calculate important properties such as CO 2 solubility for the range of expected pressure, temperature, and salinity conditions.
• Careful design of experiments combined with proxy modeling is a useful method to exploit computational resources in an efficient manner. Proxy models, when carefully built and verified, allow us to explore multitudes more scenarios compared to fullphysics simulations. Optimal solutions are arrived at much faster with proxy models, and thousands of Monte Carlo simulations can be performed within a fraction of the time taken on reservoir simulators. However, it is crucial to validate the results from proxy models as they are not constrained by physics and sometimes give results that violate fundamental physical relations, for instance, recovery factors greater than 1.

Conflicts of Interest:
The authors declare no conflict of interest.
Appendix A Table A1. Reservoir and operation parameters in the synthetic simulation model (values defined as ranges are parameters in the sensitivity analysis and uncertainty assessment).

Constraint Value Unit
Simulation time 100 years Minimum allowed producer bottom hole pressure (BHP) 2000 kPa Maximum allowed injector BHP 20,000 kPa Maximum allowed CO 2 cut in production stream 50 %