Optimization of CO2 injection using multi-scale reconstruction of composition transport

The current situation with green gas emission requires the development of low-carbon energy solutions. However, a significant part of the modern energy industry still relies on fossil fuels. To combine these two contradictory targets, we investigate a strategy based on a combination of CO2 sequestration with enhanced oil recovery (EOR) in the hydrocarbon reservoirs. In such technology, the development of miscibility is the most attractive strategy from both technological and economic aspects. Modeling of this process involves solving complex nonlinear problem describing compositional flow and transport in highly heterogeneous porous media. An accurate capture of the miscibility development usually requires an extensive number of components to be present in the compositional problem which makes simulation run-time prohibitive for optimization. Here, we apply a multi-scale reconstructing of compositional transport to the optimization of CO2 injection. In this approach, a prolongation operator, based on the parametrization of injection and production tie-lines, is constructed following the fractional flow theory. This operator is tabulated as a function of pressure and pseudo-composition which then is used in the operator-based linearization (OBL) framework for simulation. As a result, a pseudo two-component solution of the multidimensional problem will match the position of trailing and leading shocks of the original problem which helps to accurately predict phase distribution. The reconstructed multicomponent solution can be used then as an effective proxy-model mimicking the behavior of the original multicomponent system. Next, we use this proxy-model in the optimization procedure which helps to improve the performance of the process several fold. An additional benefit of the proposed methodology is based on the fact that important technological features of CO2 injection process can be captured with lower degrees of freedom which makes the optimization solution more feasible.


Introduction
Greenhouse gas emission together with a high demand of energy has long been a concern of contemporary society. Near-miscible CO 2 injection is among the most efficient strategies for a tertiary recovery of oil [16]; it can also effectively reduce carbon emissions. The produced hydrocarbons can be seen as a low-carbon fuel due to the significant amount of CO 2 left in the subsurface as the result of the EOR application. Nevertheless, the heterogeneity of subsurface with complex multi-scale characteristics requires a suitable and highly resolved model to comprehend the details of flow and interactions with the subsurface.
The current economic situation, especially low oil price and formidable cost of CO 2 , introduces extra challenges on applying a miscible gas injection. However, in combined objective of enhanced oil recovery and CO 2 sequestration, the development of miscibility may become the most attractive strategy from both technological and economic points of view. In addition, effective miscible injection can increase the storage capacity for CO 2 sequestration in virgin or depleted hydrocarbon fields. It is quite important to develop a plausible techno-economic model to meet the combined goals of oil recovery and carbon dioxide sequestration. This serves as a primary motivation for this study.
To simulate the miscible gas injection process, compositional modeling is inevitably employed. Compositional models require numerical solution of nonlinear equations that involve mass conservation of different components and thermodynamic equilibrium. The phase behavior of multiphase multi-component mixtures is usually resolved by applying an equation of state (EoS) [2,3]. A nearmiscible gas injection process usually involves a large number of species in solution, which significantly degrades simulation performance. In addition, in nonlinear iterations, thermodynamic equilibrium should be enforced in every grid block to check the phase behavior of the mixture; this adds to the performance penalty [10].
Thermodynamic equilibrium usually consists of two stages: a phase stability test [20] and flash calculation [21]. Various EoS have been used to represent thermodynamic equilibrium in a hydrocarbon mixture, starting with the classic cubic EoS [26,31]. However, the growing accuracy of reservoir fluid characterization and better recognition of complex physical processes involving component interactions requires an application of a more complicated EoS, such as statistical associating fluid theory (SAFT) [1] or cubic-plus-association (CPA) [14]. In addition, coupling with chemical reactions requires a combination of thermodynamic and chemical equilibria [17,24]. This can significantly increase the cost of phase-behavior computations in compositional simulation [34].
Several efforts have been made to improve the performance of compositional reservoir simulators by improving phase-behavior computations [9,23,28,35], spatial coarsening of compositional models [8,30], or reformulation of compositional nonlinear problem [38]. In this work, a newly proposed multi-scale reconstruction in physics (MSRP) approach by [6] is utilized for production optimization. The algebraic multi-scale (AMS) approach was initially proposed to solve an elliptic flow problem in [11]. Several extensions of this method have been successfully developed. However, most of the AMS methods were focused exclusively on the flow solver and did not address the transport problem, except [39], where an adaptive multi-scale finite volume method was proposed to accelerate the transport solver. On the basis of these ideas, an MSRP method for reconstruction of the compositional transport problem with an arbitrary number of components was developed in [6].
This approach suggests a two-stage reconstruction, where, at the first stage, the boundary of a two-phase region is recovered, and the detailed solution in the two-phase region is reconstructed in the second stage. The MSRP approach utilized an operator-based linearization (OBL) technique proposed in [36]. In the OBL method, the terms of the discretized governing equations are factorized into space-and state-dependent operators. The state-dependent operators are adaptively discretized in the parameter space of the problem, and multi-linear interpolation is applied for continuous representation [13]. This formulation helps to avoid the performance issues associated with an accurate phase-split evaluation and reduces the nonlinearity of the problem. Recently, this approach was extended for adaptive parametrization of thermal-compositional problems with buoyancy [12].
The original study of the MSRP method was limited to isothermal two-phase flow with fixed phase-equilibrium ratios (K values) [5]. In this work, we introduce an application of MSRP using the Peng-Robinson equation of state [26]. Due to the strong nonlinearity of the CO 2 injection system, constrained nonlinear optimization strategy is utilized to determine the optimal production scenario. For production optimization, we used only the first-stage MSRP reconstruction as a physics-based proxy model and compare its result with optimization of the full compositional solution. Both approaches were compared using an idealized conceptual model with growing optimization complexity.

Model description
In this section, a concise simulation framework based on [36] is presented.

Compositional framework
For simplicity, the thermal changes, capillarity, gravity, and diffusion are neglected in the following description. The general mass conservation equation for component i in the two-phase compositional problem is defined as follows: In (1), t is time, φ is the porosity of the reservoir, ρ j is molar phase density, S j is phase saturation, x i,j is the mole fraction of component i in phase j , q j is the source or sink term of phase j , and N c is number of the components. The Darcy velocity u j is defined as follows: where K is absolute permeability, k rj is the relative permeability of phase j , μ j is viscosity of phase j , and p is pressure. The equilibrium relations between oil and gas phase are required to close the system as follows: wheref i,o andf i,g are the fugacities for the component i in oil phase and gas phase, respectively. Fugacity is a function of pressure (p), temperature (T ), and phase compositions (x i,j ), which are determined by EoS-based flash computations. Additional equations are given as follows to close the system of governing equations: The overall composition of i component can be expressed as follows: where v j is the molar fraction of the phase j (o, g). Solving the Eq. (3) is a procedure called multi-phase flash [21], which will provides phase composition x i,j and phase fraction v j . Finally, the phase saturation s j can be found from the following: Applying two-point finite volume in space and backward Euler in time discretizations, the general mass conservation equation is written as follows: where V is total control volume and L represents the interface which connects the control volume with another grid blocks. In the simplified assumptions mentioned above, φ is porosity, l becomes a pressure difference between two connected grid blocks. Finally, T l j is the transmissibility of phase j . In our conceptual model, we ignored gravity, capillarity, and thermal variations focusing mostly on compositional effects. These assumptions were applied to simplify the analysis of optimization results which focus mostly on compositional effects in CO 2 injection process for EOR and sequestration. All these phenomena may increase the complexity of proxy-model, which will be considered in our future research. Notice that the OBL approach was already extended for problems with buoyancy in [12].

Operator-Based Linearization
The multi-scale technique is implemented on the basis of an OBL approach proposed by [36]. To apply OBL, the discretized mass conservation equation (8) is written in the following residual form: The operators in Eq. (9) are defined as follows: In Eqs. (10) to (14), c r is rock compressibility and T ab is the transmissibility between grid blocks. The vector u contains well-control variables, ω is the set of state variables and ξ are the set of spatial coordinates. In addition, α i is the accumulation operator, β i is the flux operator, and θ i is the source/sink operator. The OBL approach is based on a simplified representation of the nonlinear operators in the parameter space of the simulation problem. For an isothermal reservoir simulation, the parameter space is defined by the range of pressure (p) between injection and production conditions and overall compositional (z i ) range from 0 to 1. The fully implicit method (FIM) is utilized for time approximation, and Newton-Raphson method is applied to solve the governing equation Eq. (9) based on the set of nonlinear unknowns.

Multi-scale compositional transport
A solution of a compositional transport problem can be shown in a phase diagram by the solution path in compositional space, which defines the compositional changes between the initial and injection mixtures. Conservation principles and fractional flow theory form the foundation for the general solution method [22]. The compositional path of the conventional compositional problem for gas injection process, when the injection mixture is a single-phase gas and initial fluid is a single-phase liquid, always results in two shocks (leading and trailing shocks) between singleand two-phase regions. In a ternary diagram (Fig. 1a), it is presented as yellow lines connecting the initial oil and injected gas composition.
The shocks between single-and two-phase regions are always aligned along two key tie-lines (black-dashed lines) defined by liquid x i and vapor y i fractions of each  15)). Figure 1b shows the projection of the solution to fractional flow curve for CO 2 component with leading and trailing shocks (yellow) connecting tangent points on initial (red) and injection (blue) fractional flow curves respectively [22]. Note that these curves corresponds with the injection and initial tie-lines in Fig. 1a following the relation: The proposed multi-scale compositional transport approach consists of two stages [5]. The first stage utilizes the set of restriction-prolongation operators for reconstructing twophase boundaries (the trailing and leading shocks). The restriction here reduces the n c − 1 transport equations to a single equation with a special flux operator based on the pseudo-fractional flow curve. In the second stage, the set of restriction-prolongation operators is applied in the two-phase region to reconstruct the solution structure of the two-phase displacement. This stage is based on the invariance of two-phase solutions in tie-line space reported in [33] and adapted for practice in [35]. The proxy model for compositional simulation, utilized in this work, uses the first-stage multi-scale reconstruction from [5]. A restriction operator combines two fractionalflow curves for injection and production tie-lines (red and blue curve from Fig. 1b), defined as follows: The equivalent fractional flow curve (green in Fig. 2), serving as the restriction operator, is constructed by taking a convex hull on the union of both curves.
This means that the green line in Fig constructed pseudo-fractional flow curve. Notice that by structure, this system is very close to the conventional binary compositional problem. Once the solution of the restricted system is found, the full system is reconstructed based on the prolongation operator. This operator applies interpolation for all components in the solution between initial and injection compositions using the solution of the restricted system κ(z R ) (corresponds to the CO 2 component in this example) as an indicator: Here, κ is the interpolation-prolongation operator, z R is the restricted solution, and I is the piecewise linear interpolation function. Referring to this linear interpolation, the transport solution of other components in the multicomponent system is reconstructed and used as a proxy model in place of the full compositional model. Notice that this system can accurately predict only the boundaries of the twophase region and their dynamic propagation in space; for a really accurate solution, the second-stage multi-scale reconstruction should be applied [5].

Economic model
The techno-economic model is applied to evaluate the economics of a combined CO 2 EOR and sequestration application. Several economic studies of CO 2 injection processes have been performed in [4,15,29,32]. McCoy and Rubin [19] proposed several regression equations for assessment of the capital cost of CO 2 injection projects, which are validated in [37] and [4]. Referring to [32], this techno-economic model uses simulation input data and oil production rate, gas injection rate, and bottom hole pressure (BHP) to define different costs and revenues of the project. On the basis of reservoir-simulation data, an economic model is developed to estimate the profitability of CO 2 injection for enhanced oil recovery (EOR) and CO 2 sequestration, which will reflect on the net present value (NPV). The general economic parameters of a CO 2 injection process are listed in Fig. 3.
This figure shows that the cost of a CO 2 injection project can be divided into two parts, which are capital cost and operational cost. Dominant revenues from the gas injection project mainly originate from oil sales and carbon sequestration incentives. A previous economic study of CO 2 injection projects [15] indicates that CO 2 purchasing cost is one of the most sensitive parameters when NPV is evaluated.
In this work, we identify that CO 2 processing cost has a similar impact on NPV as CO 2 purchasing cost. The CO 2 processing cost model in this work is based on [32], and is expressed in terms of the pump capital cost as follows: where W p is pumping power requirement, which is expressed in kW, which in turn varies with CO 2 injection pressure. Other parameters in the economic model are listed in Table 1. Some of them are obtained by introducing the regression equations listed in [19], such as those for well engineering cost and CO 2 processing equipment cost.

Numerical results
In this section, we demonstrate the comparison between solutions of the proxy model and the full compositional model. Here, we limit our investigation to a conceptual 1D reservoir model for simplicity of the optimization results interpretation. In this model, the injection well on the left operates at a constant gas rate when the production well is controlled by BHP which serves as a control variable for optimization. Figure 4 shows the estricted solution z R , which yields the shock reconstruction curves for simulation results for the growing BHP at the production well. All simulation results are shown for the model with parameters specified in Appendix A. The K-value table in this work is obtained from the embedded constant composition expansion (CCE) experiments in [7] based on the PR EoS, which is shown in Tab 10. It is clear that the K value system does not develop miscibility even when BHP provides the pressure at the displacement front close to the first-contact minimal miscibility pressure (FC MMP) for this system (around 126 bars at T = 373K). This happens due to the inability of the K-value model to predict miscibility accurately, since compositional dependency is not captured in this model. It can be overcome by either extension of the K value parameterization with additional degrees of freedom, e.g. [27], or incorporation of EoS-based phase behavior [2]. However, it is clear that the two-phase boundaries can be accurately represented by the restricted model for K-valuebased physics. In addition, the complexity and structure of the restricted solution are invariant with respect to the number of components present and only depends on initial and injection tie-lines in the multicomponent system (see, [5] for details). Next, the results of the restricted solution for the compositional problem based on the EoS is shown. The structure of the compositional transport solution depends on key tie-lines [22]. For the restricted solution, we follow the same strategy as before and construct the restriction operator based on combined fractional flow (16) according to the first stage of MSRP approach [5]. The solution of the restricted transport equation reconstructs the boundaries of the two-phase region using one transport equation instead of n c − 1 equations in the conventional compositional model.

Restricted solution
The results of quaternary system reconstruction are shown in Fig. 5. Here, you can see that for a high BHP value, the structure of the solution is much closer to miscibility (leading and trailing shocks stays closer to each other) than in the K-value approximation. This happens because the EoS-based phase behavior correctly represents the compositional dependency of the solution. Similar to the K-value model, the restriction stage requires the solution of only one equation instead of n c − 1, where n c is the number of components.
Next, we present the simulation results for more realistic multicomponent mixture. Here, we used the eight-component system from [38] with compositional parameters shown in Table 9 (see Appendix). In cooperating with the K values generated using [7], and shown in Table 11, the restricted solution based on the first-stage reconstruction of MSRP approach, is present in Fig. 6.

Prolongation of proxy model
Here, we illustrate the construction of the proxy model using an interpolation-based prolongation operator Eq.  (Tables 9 and 11). This proxy model still cannot recover the full solution of the eight-component model and the secondstage MSRP approach should be employed to reconstruct the solution in two-phase region (see details in [5]).
The prolongation yields a full compositional solution at every grid block, which then can be used in a multiphase flash procedure to predict phase behavior. This provides an accurate reconstruction of the two-phase boundaries. In our proxy model, we are using this prediction to compute phase rates at wells. As a result, this proxy model will be applied to evaluate economic performance of the CO 2 injection project together with full eight-component model in terms of NPV. Figure 11 shows the transport solution for both fully compositional and proxy models for different BHP controls at the production well and a fixed rate at the gas injection well. It is clear that with increasing pressure, both models capture the development of miscibility since the BHP control at the production rate will apparently control the pressure at the displacement front. Due to the development of miscibility, the leading and trailing shocks get closer to each other and the displacement efficiency grows. Next, we investigate optimal production strategies for this model. In the optimization stage, the full four-component system together with the proxy two-component system is used to determine oil production. Net present value (NPV) is used as an indicator to estimate the economic profitability of the project. The simulation time is divided into several periods where changes in BHP at the production well are applied. Here, we make sure that the time period for simulation covers the breakthrough of the trailing shock at the lower limit of pressure. Next, we estimate the optimal production strategy with a different numbers of control variables.

NPV with a limited number of control parameters
The NPV distribution as a function of a single BHP control is evaluated here. We compare the NPV curve vs. BHP control for both the proxy and the full compositional model. The simulation time is defined to be long enough for the breakthrough of both leading and trailing shocks of the solution. The NPV plotted as a function of control BHP is shown in Fig. 12.
Here, the green solid curve is the NPV results from the full four-component model, and the red-dashed curve is the NPV results from the proxy model. While there are some discrepancies in the proxy solution due to the limited application of the MSRP (only first stage of reconstruction), To reduce the differences in NPV evaluation, the second stage of the MSRP reconstruction can always be performed. Notice, that this behavior is expected. In most of current gas injection projects, the cost of CO 2 remains a major factor for project economy. In our model, the simplified physical assumptions and CO 2 sequestration credits yield the development Fig. 12 NPV with one control parameter of miscibility as an optimal strategy to improve NPV. The BHP corresponding to the maximum value of NPV in Fig. 12 is located close to the minimal miscible pressure of this system which is followed by the NPV reduction due to the growing cost of pumping. The majority of proxy models for gas injection problems poorly approximate the development of miscibility which requires an application of high-fidelity compositional models. We demonstrate that our proxy model is able to capture the near-miscible behavior and correctly identify the maximum of the NPV.
Next, we introduce two simulation time periods and two control variables (BHP 1 and BHP 2 ) for NPV evaluation. For this situation, the prediction of optimal controls for both periods is not as simple as in the previous example. Performing an exhaustive search in the space of control variables, we evaluate the NPV function (shown in Fig. 13). While the NPV function is different for the proxy and the full model, the maximum NPV is reached at the same control values, i.e., around BHP 1 = 95 bars and BHP 2 = 118 bars. These values are conditioned by the obvious strategy for production controls when in the first-time interval, the lower BHP at the production well provides the near-miscible pressure at the displacement front. In the second-time interval, the higher pressure at the production well provides nearmiscible pressure until the breakthrough to the production well. The near-miscible strategy is optimal since it maximizes the oil recovery and sequestrationof CO 2 . In order to illustrate optimization procedure, the optimization trajectory is constructed. Figure 14 gives the result of the optimization trajectory for the full physics model and proxy model respectively. For the proxy model, the number of optimization steps is less than that of the full physics model. In addition, it is faster for the proxy model to acquire a near optimal NPV result with the same optimization strategy. Table 2 presents more details of both optimization runs.
Next, we look into the form of the NPV function for the more realistic eight-component system with compositions and corresponds thermodynamic characteristics shown in the Appendix (Tables 7 and 9). The objective function for two control variables is shown in Fig. 15.
The pressure interval is corresponding to the lower and upper limit of the BHP values of the production wells, which range from 80 to 160 bars. In this eight-component system, given the similar optimal pressure sets for both proxy and full physics model, both models capture the similar highest NPV value. The result is shown in the Table 3.

Optimization with multiple controls
Next, we apply production optimization based on five control variables (BHPs) corresponding to five time periods in the simulation. In this study, we use the "fmincon" function from the Matlab optimization toolbox [18]. In "fmincon," the "sqp" algorithm has been chosen. The optimizer is utilized to provide BHP controls at each time period and obtain an optimal NPV result during the CO 2 injection process. All BHP controls were bounded by BHP min = 60 bars and BHP max = 140 bars. Note that the expected optimal strategy should include a gradual increase of BHP at each consecutive control interval to provide near-miscible conditions at the displacement front. We test several initial guesses for the optimization with five control parameters. For this number of controls, several local minima can exist and the optimizer struggles with finding a single global extremum. However, based on the  Table 4. You can see that the proxy model performed fewer iterations and obtained a similar NPV.
In addition, we perform two more optimization runs with different initial guesses when all BHP controls have been set to 70 bars and 100 bars respectively. The results can also be seen in Table 4. In these optimization runs, both models cannot converge to the same optimal strategy but get close to it. The proxy model performs quite robustly and proves to be applicable for optimization of gas injection process in the idealistic reservoir.

Conclusions
In this work, we extend the multi-scale reconstruction in physics (MSRP) approach for the EoS-based gas injection problems. In particular, we parametrize the restriction operator of the first-stage MSRP reconstruction in the pressure interval and obtain the restricted solution using the operator-based linearization framework. The restricted solution was prolongated to the full compositional solution using interpolation operator. The obtained proxy model can accurately predict the boundaries of the two-phase region and has been utilized in this work for production optimization in a simplified physical assumptions of the forward problem. Referring to previous economic assessments of CO 2 injection projects, a techno-economic model has been developed to analyze the revenues of CO 2 injection for the combined objective of EOR and sequestration. A general application of the proposed proxy model for optimization of gas injection process is demonstrated in this study. Starting with a limited number of controls, we show that the objective function of the full physics compositional model and the proposed proxy model share similar extrema for a limited number of control parameters. To test the robustness of the proposed proxy model in relatively complicated cases, the general form of the objective function was evaluated for a limited number of control parameters. Based on these evaluations, we demonstrate that both full physics and proxy models share similar extrema.
In addition, a constrained nonlinear optimization is applied to determine an optimal production strategy for the gas injection operation in the simplified physical setting. For an optimization model with more control parameters, when the initial guess of controls is near the optimal solution, we show that both full physics model and proxy model converge to similar optimal solution. For arbitrary initial guesses, the converged optimal strategy may differ between the proxy and the full compositional models due to a local extrema of both objective functions.
Through the optimization process in four-component system, we have shown that by providing the optimizer with the same input parameters in both the full physics model and the proxy model, the optimal solution with the proposed proxy model is usually more feasible (takes less iterations) than the full physics model. In addition, the forward simulation of the proposed proxy model is significantly cheaper (proportional to the reduction in the number of components) than the full-physics model and becomes comparable with the conventional black oil model. In our future work, we will extend the proposed model for a more realistic situation involving more governing physics.   Tables 10 and 11 for details).
For the K-value model, we perform the CCE using PVTi module [7], where we generate a K-value table corresponding to given initial compositions in Tables 8 and 9. The K-value table is present as a function of pressure with three pressure values employed (see Table 10 and 11 for details).