Optimization of engineered water injection performance in heterogeneous carbonates: a numerical study on a sector model

One of the emerging technologies for boosting oil recovery in both sandstone and carbonate reservoirs is engineered/low-salinity water injection (EWI/LSWI). In this paper, optimization of engineered water injection is investigated using three synthetic sector models representing homogeneous, heterogeneous with channeling, and heterogeneous with gravity underride reservoirs. Both oil recovery and net present value were investigated as objective functions for the study. Eighteen design parameters were selected for the study including reservoir, operational, and economic parameters. Response Surface Methodology and Designed Exploration and Controlled Evolution algorithms were implemented for sensitivity analysis and optimization studies, respectively. The study highlighted that NPV is more representative as an objective function compared to oil recovery where the three optimized models have about similar oil recovery, but different NPVs. The sensitivity analysis showed that oil price, tax rate, and initial oil saturation are the three most influential design parameters on the net present value for the three models investigated. Moreover, the findings showed that developing the gravity underride model requires more attention as being the most sensitive model with 13 influential design parameters. The optimization study highlighted that secondary EWI is recommended to achieve the best profitability out of the three models. However, a high maximum exposure is expected due to the capital and operational costs related to early EWI application. This study is one of the very few that discusses the economic aspect of EWI while incorporating the complexity of geochemical reactions and the heterogeneity of carbonates.


Introduction
Significant improvements in oil recovery from carbonates have been observed through adjusting the ionic composition of the injected brine (Webb et al. 2005;Zhang et al. 2007;Jerauld et al. 2008;Austad et al. 2008;Yousef et al. 2011; 1 3 Gupta and Mohanty 2011;Chandrasekhar and Mohanty 2013;Almeida et al. 2018; Ajayi and Gupta 2019;. While the mechanism underlying this observed increase in oil recovery is still unclear, most of the researchers have consensus on wettability alteration as the main reason in carbonates (Strand et al. 2006;Zhang et al. 2007; Austad et al. 2008;Yousef et al. 2011;Chandrasekhar and Mohanty 2013). In the context of this work, the process of diluting the injected brine while maintaining the same proportions of different ions is termed as low-salinity water injection (LSWI). On the other hand, the process of modifying the ionic composition, by either increasing the concentration of the ions (hardening) or decreasing the concentration of the ions (softening), is defined as engineered water injection (EWI).
Few numerical works in the literature capture the wettability alteration observed in the laboratory by LSWI/ EWI, especially in carbonate reservoirs. This is due to the complexity of the geochemical reactions involved as well as the high heterogeneity of these reservoirs. For field-scale predictions of oil recovery by LSWI, Jerauld et al. (2008) published a numerical model that is based on interpolation between two sets of relative permeability and capillary pressure curves. Qiao et al. (2016) proposed a mechanistic model of LSWI through incorporating the geochemical interactions between crude oil/brine/rock (COBR) using the PennSim Reservoir Simulator. Their results were in agreement with the effluent-ions concentration measured from the chromatographic test conducted using limestone cores. Ghedan et al. (2017) used the concentration of sulfate to model irreversible wettability, which they claimed is applicable to both sandstone and carbonate reservoirs. Awolayo et al. (2017) used both Ca 2+ and Mg 2+ on the rock surface charge as interpolants to history match published coreflood data. Yutkin et al. (2018) presented a simple ion-complexing model to predict the calcite-surface charge as a function of diluted brine composition that is saturated with carbon dioxide. Adegbite et al. (2018a) proposed a geochemical model for engineered water injection in carbonate reservoirs, which was validated through history matching of recent coreflood data. This proposed model will be thoroughly described later in this paper.
Furthermore, few studies incorporated carbonate heterogeneity through examining the effect of channeling and gravity underride on oil recovery. Cho et al. (2015) observed that during a water-alternating-gas (WAG) injection process, the injected water settled at the bottom of the model while the gas rose to the top. Also, Chang et al. (1994) observed that the flow pattern was channeling when the permeability fields had a large Dykstra Parson's coefficient (V DP ) and a large correlation length. Adegbite et al. (2018b) investigated the effect of EWI on oil recovery from carbonate reservoir with channeling and gravity underride problems. The results highlighted the adverse effect of heterogeneity on oil recovery by EWI as opposed to a homogeneous model. The previously described complexity of geochemical reactions and heterogeneity of carbonates necessitate further sensitivity and optimization studies on EWI in carbonates. P'azman (1986) defined Experimental Design as a set of observations that are implemented to achieve a certain objective function. One or more design variables are changed to observe the effect on a response variable. Ghorbani (2008) showed that these design variables are usually divided into state parameters involving uncertain parameters that cannot be controlled and decision parameters that can be controlled. Mason et al. (1989) defined the Response Surface Methodology (RSM) as a geometric demonstration of the objective function. The idea behind these two methods, Experimental Design and Response Surface Methodology, is to build up the response surface that truly represents the response of the objective function. Damsleth et al. (1992) and Dejean and Blanc (1999) suggested different approaches involving Response Surface Methodology and Monte Carlo simulation. Another method is the Designed Exploration and Controlled Evolution (DECE) optimization method, which reservoir engineers commonly use to solve history matching or optimization problems. DECE can be described as an iterative optimization process that first applies a designed exploration stage followed by a controlled evolution stage using a specific objective function (Nguyen et al. 2016).
Very few studies in the literature discussed optimization of LSWI/EWI. Al-Shalabi et al. (2014) investigated the optimization of low-salinity water injection in carbonate reservoir using the cumulative oil recovery as the objective function. They used seven uncertain and decision design variables by means of a 5-spot LSWI pilot models using the UTCHEM simulator. They achieved an optimum LSWI design of 73.57% cumulative oil recovery for the design parameters tested. Nguyen et al. (2016) investigated the field-scale optimization of CO 2 -LSWAG under geological uncertainties using clay distribution and Designed Exploration and Controlled Evolution (DECE) optimization method. According to them, DECE will efficiently explore the solution space and find the location of the injector/producer pattern that maximizes the ultimate oil recovery. They achieved a recovery of 43.5% for the optimal solution compared to 30.2% for the base case. More information about the EWI/ LSWI techniques including experimental, numerical, and field works can be found elsewhere (Al-Shalabi et al 2015a, b; Al-Shalabi and Sepehrnoori 2017).
As seen from the above-mentioned optimization studies, the main focus was on the technical aspect of maximizing oil recovery rather than the economic aspect. Profitability is the major driver and the design governor of any EOR project (Hite et al. 2005). Many projects were deemed profitable in the oil and gas industry for many years because of the increase in oil demand as well as prices, which enabled companies to invest in new technology and operations. However, the drop in oil prices to about $ 20 per barrel has caused difficulties in guaranteeing the profitability of any project (International Energy Agency 2019). Bondor (1993) reported that EOR projects should be subjected to strict criteria than previous practices. For a project to be successful, the economic study is performed together with engineering design and performances. The latter ensures that extra attention is given to factors that have great impact on profitability to keep the project within budget (Hite et al. 2005).
Several economic indicators are used to measure profitability and to rank multiple projects. These indicators include net present value (NPV), the time value of money, differential cash flow, internal rate of return (IRR), and payback period. In this work, the NPV is mainly targeted through optimization with some insight into payback period and maximum exposure. The NPV is used to make key decisions when cash flows are received at different times. It is defined as the total present value (PV) of time series of cash flow and is expressed as follows: and present value is given by: where t is the time of the cash flow, I is the discount rate (interest rate), and R t is the net cash flow (positive for inflow and negative for outflow) at time t. A project is said to be profitable if the NPV is greater than zero and unprofitable if it is less than zero (Albright and Winston 2012). The payback period is another economic indicator that describes the duration at which the initial cash outflow of an investment is expected to be recovered from the cash inflows govern by the investment, and it is expressed as follows (De Jonge 2016): Maximum exposure can also be used to describe the project profitability, and it is defined as the maximum negative cash flow on a project (Norwegian Petroleum Directorate 2004).
This study is one of the very few that investigates numerical optimization of engineered water injection from an economic prospective while incorporating the complexity of geochemical reactions and the heterogeneity of carbonates. This optimization incorporates the complexity of geochemical reactions through using a mechanistic geochemical model. Also, the heterogeneity of carbonates The cost of Project∕Investment Annual Cash Inflows .
is incorporated through capturing channeling and gravity underride scenarios. In addition, the objective function used for the optimization study was selected based on technical aspect (maximizing oil recovery) as well as economic aspect (maximizing net present value). The paper includes descriptions of experimental input data, setup of simulation models and simulation data, engineered water injection model, and the EWI sensitivity and optimization approach.

Model setup and simulation data
The CMG-GEM™ compositional simulator was used in this study for setting up sector models to predict FW and EW injections in carbonate reservoirs. The data used in these sector models are based on the experimental work of Chandrasekhar and Mohanty (2013). Several corefloods were conducted by Chandrasekhar and Mohanty (2013) among which two corefloods were utilized for the study, formation water (FW) and engineered water (EW) using magnesium and calcium with four times spiked sulfate (MgCa4S). The latter two corefloods were history matched in our previous works (Adegbite et al. 2017(Adegbite et al. , 2018a. Further information about the laboratory studies can be found in the work of Chandrasekhar and Mohanty (2013). Also, more details about the simulation and history matching of these laboratory studies can be found elsewhere (Adegbite 2018;Adegbite et al. 2017 and2018a). Three synthetic 5-spot sector models were created by upscaling the history matched corefloods of Chandrasekhar and Mohanty (2013). The upscaling of corefloods to sector models was performed through maintaining the same flow functions obtained at the core scale and applying them to the sector scale. These flow functions include relative permeability and capillary pressure curves. For the static properties including porosity and permeability, the average value at the core scale was maintained and then, was distributed at the sector scale. The heterogeneous carbonate core plugs had an average porosity of 19.6% and an average brine absolute permeability of 40.6 mD. The sequential Gaussian simulation (SGS) algorithm in the Petrel software was utilized for creating the permeability and porosity distribution maps using their average values (Petrel Technical Manual 2017). In the SGS method, a variogram correlation and a probability distribution are needed for the simulations. A spherical variogram was utilized which represents a medium level of heterogeneity between Gaussian and exponential variograms. In addition, log-normal and normal distributions were assigned for permeability and porosity, respectively. In these probability distributions, a mean and a standard deviation are needed to define the distribution. The mean was taken as the average core value, and the standard deviation was estimated using a Dykstra Parson's coefficient (V DP ) of 0.8 for both permeability and porosity heterogeneous maps. Hence, three sector models were created including homogeneous model (base case), heterogeneous model with permeability channeling, and heterogeneous model with gravity underride (Fig. 1).
The models were set up in a 3D Cartesian grid system with 17 × 17 × 6 gridblocks. It is worth mentioning that the number of gridblocks used was based on a grid refinement analysis on each of the three sector models. This was performed by comparing the oil recovery for the 1734 gridblocks base model (Table 1) with varying total gridblocks of 500, 1000, 1575, 2000, 4000, and 5000. In addition, this analysis was performed in the tertiary mode of injection (formation water injection first, followed  by engineered water injection) to check for numerical accuracy and to capture the physics of both formation and engineered water. The results, as depicted in Fig. 2a-c, show that the difference in oil recovery for all the six models compared with the base case is less than 3%. Moreover, an additional study was performed to investigate the impact of physical dispersion on sulfate-ion concentration and transport during engineered water injection (Fig. 2d). The latter figure shows that the changes in velocity-dependent longitudinal and transverse dispersivities have no major significance on sulfate wave propagation during engineered water injection. The simulation runs were performed in both secondary and tertiary injection modes. Tables 1 and 2 describe the dimensions and inputs for the sector models. For the channeling and gravity underride models, 10 ft longitudinal dispersivity and 1 ft transverse dispersivity were used. Table 3 shows the calculations of the capillary number (8.59 × 10 −8 ) and pressure gradient applied (1.39 psi/ft) to ensure that neither the critical capillary number nor fracturing pressure is exceeded. One should note that the capillary number (N c ) was calculated using the following equation (Jin 1995): where k is the formation permeability, σ is the water-oil interfacial tension, ΔP is the pressure drop between the injector and the producer wells, and L is the distance between the injector and the producer wells. Marhoun (1998) and Ahmed (1989) correlations were used to estimate oil formation volume factor (B o ) (1.862 bbl/STB), which was later used to determine original oil in place (OOIP) ("Appendix 1"). Other simulation data were taken from Chandrasekhar and Mohanty (2013).

Engineered water injection model used
The mechanistic geochemical model used in this study to capture the complex reactions in carbonates is based on an anion exchange reaction. The latter reaction occurs on the positively charged surface sites in carbonates, which act as exchangers. Brady and Thyne (2016) supported the presence of these positive sites, which render the net positive charge usually found on carbonate rocks. Adegbite et al. (2018a) proposed that during EWI, oil carboxylic group undergoes anion exchange reaction with the sulfate-ion present in the EW. This alters the wettability of the carbonate rocks toward a more water-wet condition and hence leads to the release of oil ganglia (Fig. 3a). The governing reaction is given by:

Fig. 3
Wettability alteration in carbonate rocks. a EWI mechanism for wettability alteration proposed by Adegbite et al. (2018a). b LSWI mechanism for wettability alteration proposed by Zhang et al. (2006) where CH 3 COO − represents oil carboxylic group and X represents the exchanger. This equation was used for the geochemical modeling of oil recovery by EWI using CMG-GEM™. This proposed model differs from Zhang et al. (2007) model (Fig. 3b) where the latter necessitates the presence of other ions, namely calcium ion and/or magnesium ion, along with sulfate to cause wettability alteration. For Adegbite et al. (2018a) proposed model, SO 4 2− is the key contributor that makes the rock more water wet and thus increases the oil recovery. One should note that although the interpolating parameter is mainly expressed in terms of SO 4 2− ion, the contribution of other ions and particularly Mg 2+ and Ca 2+ ions is indirectly incorporated through affecting sulfate concentration along with the other incorporated geochemical reactions. Sulfate equivalent fraction, Σ(SO 4 -X 2 ), as the interpolating factor was used to interpolate between two sets of relative permeability and capillary pressure curves (Figs. 4 and 5). The set of curves used in Fig. 4 represents a mixed-wet rock; however, those in Fig. 5 represent a water-wet rock. It is worth mentioning that these two sets of relative permeability and capillary pressure curves are different as they represent two wettability states; a mixed-to-oil wet state (initial state) and a water-wet state (final state). The initial state of mixed-to-oil wet usually prevails as an initial condition in carbonates, while the final state of water-wet state represents the wettability alteration due to the use of engineered water injection. These both sets of curves were obtained from our previous works through history matching the experimental corefloods of Chandrasekhar and Mohanty (2013). More

Fig. 4
Relative permeability and capillary imbibition curves for FW (Adegbite et al. 2018a). a FW relative permeability curves (k r ). b FW imbibition capillary pressure curve (P c )

Fig. 5
Relative permeability and capillary imbibition curves for EW (Adegbite et al. 2018a). a EW relative permeability curves (k r ). b EW imbibition capillary pressure curve (P c ) information about these curves and the history matching process can be found elsewhere (Adegbite et al. 2018a).
The latter set of curves will be utilized for injecting engineered water in both secondary and tertiary modes of injection. In this work, different water compositions with varying sulfate concentrations will be utilized. Depending on the sulfate concentration in the engineered water injection, interpolation occurs between the mixed-wet and the waterwet sets using sulfate equivalent fraction as an interpolation factor. It is worth mentioning that the same set of curves were used across the three sector models.

EWI sensitivity and optimization approach
Response Surface Methodology (RSM) and Designed Exploration and Controlled Evolution (DECE) algorithms have been widely used in optimizing several chemical processes (Al-Shalabi 2014; Li et al. 2015). This is due to the importance of these techniques in reducing the uncertainties involved in any project. In this work, RSM and DECE were used to address the sensitivity and optimization of EWI technology in carbonate reservoirs using both oil recovery and net present value (NPV) as objective functions (response variables). Several steps were considered for the proposed sensitivity and optimization approach: 1. Determining the response variable to be tested as well as the design parameters that affect the response variable with their respective range of values (maximum and minimum).

Using Response Surface Methodology (RSM) and
Experimental Design to generate and run different simu-lation cases, and then to build a proxy (approximation) model which represents the initial reservoir model. 3. Evaluating the sensitivity of each design parameter on the response variable, after which insignificant parameters were screened out. It is worth mentioning that the latter tools calculate the contribution (%) of each design parameter on the response variable and accordingly a cutoff (%) is determined to exclude the insignificant parameters. 4. Performing numerical optimization using DECE optimizer to obtain the optimal design using NPV as objective function.
A description on the selected response variable as well as design variables is presented below.

Response variable
Two response variables were initially considered for the study as was previously mentioned; however, later only NPV was selected as a main parameter for optimization. Maximizing NPV was the main objective function used. Details on eliminating oil recovery as a response variable will be provided later in the results and discussion section. The economic input data used in NPV calculations for the sensitivity and optimization studies are shown in Table 4.
These economic input data are divided into three categories including capital costs, operating costs, and tax and economic. Capital costs involve facilities and equipment as well as intangible costs. Operating costs include FW and EW injection operating costs, produced water, oil treatment, and overhead costs. Tax and economic include royalty, inflation rate, oil price escalation rate, operating cost escalation rate, tax rate, discount rate, and oil price. It is important to note that these economic data, especially the royalty and tax rate, are applicable to the Middle East region (Ahmed 2018). Also, the oil price of $61 per barrel was chosen as a base case that represents the oil price at the time of conducting this study. A more detailed justification of the chosen values is presented below in the Justification of ranges used for design parameters section.

Design parameters
Eighteen design parameters were selected for the study as listed in Table 5. Different aspects were considered in the selection of these 18 design parameters including the technique (engineered water injection), the reservoir at which the technique is applied (carbonates), and the objective function selected (maximizing NPV). Therefore, the design parameters are divided into reservoir, operational, and economic parameters. The reservoir parameters are EW residual oil saturation, FW residual oil saturation, selectivity coefficient, irreducible water saturation, anion exchange capacity (AEC), calcite volume fraction, dolomite volume fraction, transverse and longitudinal dispersivities, cross-flow (K v /K h ratio), and reservoir heterogeneity (V DP ). The latter parameters are naturally existing and have a potential of affecting the performance of engineered water injection. Operational parameters, which usually the field operator has a control over, include injected water pore volumes of EW and FW as well as water composition with varying sulfate concentration for EW and FW. Economic parameters, which are controlled by the oil and gas market and other constraints, include oil price, tax rate, and royalty.

Justification of ranges used for design parameters
It is worth mentioning that the range of values selected for each design parameter was based on typical values from the literature and our expertise that reflect the nature of the targeted carbonate reservoirs and the applied engineered water injection technique. The range of values for the tax and royalty rates was chosen based on the Middle East countries surveyed in the work of Ahmed (2018) (2016). However, one of the most important parameters is the desalination plant for engineered water that represent 2% of the capital cost similar to BP's percentage cost for a similar project (2b1stconsulting, 2019). This is taken together with the other cost for EW operating cost. The other economic parameters were obtained from the work of Ghorbani (2008) and Layti (2017). It is important to note that some of the selected costs reflect the 5-spot sector models used. For the non-economic parameters, FW and EW sulfate concentration ranges were chosen based on the minimum and maximum concentrations at which there was no effect on oil recovery as was observed in the laboratory work of Chandrasekhar and Mohanty (2013). Also, FW and EW residual oil saturations were chosen from the work of Chandrasekhar and Mohanty (2013). The irreducible water saturation was chosen based on the minimum and maximum values that

Results and discussion
This section presents the results and discussion of EWI numerical optimization in both secondary and tertiary injection modes for the three sector models. The workflow starts with performing a sensitivity analysis on the design parameters followed by an optimization analysis study. These two steps are thoroughly discussed below for the three sector models starting with heterogeneous with channeling, followed by heterogeneous with gravity underride, and then homogeneous models. Both NPV and oil recovery were initially considered as objective functions; however, later the NPV objective function was selected for the rest of the analysis as justified from the discussion.

Sensitivity analysis (SA)
This section investigates the sensitivity of design parameters on the three sector models by performing series of reservoir simulations using both Design of Experiment (DOE) and Response Surface Methodology (RSM) algorithms. Twolevel fractional factorial design as a type of DOE was considered at first to reduce the number of simulation runs and provide information about the relative effects of each design parameter on the response variable. Then, the RSM was used to explore the relationship between the design parameters and the objective function by running design experiment sufficiently enough to build a proxy model. The RSM was utilized to build proxy models of high-degree polynomial through performing regressions and minimizing the sum of squares error between the true response variable and the predicted response variable. The RSM utilizes different types of polynomial functions to fit through the data including intercept, linear terms, quadratic terms, and even, cubic terms. It is usually sufficient to use up to the quadratic terms. The higher the degree of the polynomial, the closer the Taylor series can approximate the true value of the response variable. An example of the full quadratic polynomial model is defined as (Myers et al. 2009): where x i and x j are independent input parameters, y is a response variable, and β0, β i , β ij , and β ii are the regression coefficients that are selected to minimize the summation of the squares of the error (ε). The proxy models created for this study are discussed later in this section and given in "Appendix 2." These proxy models were then used to determine the sensitivity of each parameter through tornado plots as well as for optimization studies.

Heterogeneous model with channeling
Two SA studies were considered for this model using two different objective functions; net present value (NPV) and oil recovery. The first objective function is the net present value that has to be maximized based on the 18 design parameters (Table 5). Likewise, the second objective function (oil recovery) has to be maximized using the same design parameters used for the NPV; however, the economic input parameters (tax rate, royalty, and oil price) were not considered. About 90 simulation runs were performed for each objective function where the range of each design variable was randomly selected using the Experiment Design algorithm.

NPV objective function
For the NPV objective function, the Response Surface Methodology (RSM) was used to explore the relationship between the design parameters and the objective function by running design experiment sufficiently enough to build a proxy model. The latter model was then used to determine the sensitivity of each parameter through a tornado plot. Figures 6 and 7 show the oil recovery and the corresponding NPV obtained for each of the simulation runs performed, respectively. One should note that NPV was calculated using Eq. (1) for each of the design experiments. It is also worth mentioning that the black color in these figures represents the base case run with parameters defined in Tables 2  and 4 for reservoir and economic inputs, respectively. Figure 6 shows that some of the simulated runs are lower in oil recovery compared to the base case; however, the majority of those runs resulted in higher recovery. Also, the latter figure depicts that some of these simulation runs were performed in the secondary mode of injection with zero injected pore volumes of formation water. On the other hand, the base-case simulation resulted in the lowest NPV as depicted in Fig. 7. Therefore, there is a potential to improve NPV within the ranges used for the different design parameters. After the simulation runs, a proxy model was generated to investigate the sensitivity of the design parameters. The proxy model generated for the NPV is given by Eq. (1B) in "Appendix 2." This equation shows the 18 design parameters considered in this study. The verification plot of the generated proxy model is shown in Fig. 8. This figure compares the NPV calculated through simulation and that predicted using the proxy model. Each simulation run is presented with a point, and the distance of each point from the 45-degree line represents the error. The figure shows that the proxy model has fair predictions of the NPV since most of the simulations are close to the 45-degree line. Therefore, the proxy model for the channeling model is reliable for further sensitivity analysis and optimization.  The results of the sensitivity analysis conducted on the proxy model for channeling is presented in a tornado chart (Fig. 9). This chart indicates the degree to which each design parameter affects the NPV. This figure enables the quantification of the relative importance of each parameter. The SA results show that oil price, tax rate, and irreducible water saturation are the three most significant parameters that affect the profitability of the EWI project. Other significant parameters are FW injected pore volume, EW residual oil saturation, royalty, FW residual oil saturation, reservoir heterogeneity (V DP ), EW sulfate concentration, and EW injected pore volume. Therefore, 10 out of the 18 design parameters are significant and will be used later in the optimization study, while the rest of parameters were screened out. It is interesting to note that the economic parameters have the highest effect on the heterogeneity model with channeling as opposed to other reservoir and operational parameters. It is important to note that in this case, a design parameter is considered significant if its effect on the objective function (NPV) is greater than 0.5%. The screening criteria were chosen based on the analysis of the statistical parameters.
Oil recovery objective function The same previously described approach was followed for the oil recovery as an objective function. Nevertheless, it is important to note that the economic parameters including oil price, tax rate, and royalty were not considered in this case. The proxy model generated for analyzing the sensitivity of oil recovery for the heterogeneous model with channeling is given by Eq. (2B) in "Appendix 2." As was previously mentioned, this equation reflects the 15 design parameters considered excluding the economic parameters.
The tornado chart for this case is presented in Fig. 10. One can note that the three most influential parameters are EWI pore volume, FWI pore volume, and EW residual oil saturation. One should note that pressure maintenance is not the only reason for EW and FW pore volumes injected being the most significant parameters, but water salinity is also important. This is supported by a higher oil recovery by engineered water (EW) as opposed to formation water (FW) for the same pore volumes injected (Adegbite et al. 2018a). Actually, in most of the optimized scenarios, as discussed later in the optimization section, FW has lower PVs injected compared to EW. The latter address the proper injection time of EW and encourages applying EW in the secondary mode of injection. Other significant parameters are highlighted in the figure showing that 13 out of the 15 design parameters tested are significant.
It is interesting to highlight the differences between oil recovery and NPV as two objective functions. In oil recovery objective function, operational and reservoir parameters become important, e.g., EW injected pore volume and EW residual oil saturation. On the other hand, for NPV as an objective function, other economic parameters turn out to be significant, e.g., oil price and tax rate. In addition, the proxy model generated using the NPV objective function represents a first-order multivariable polynomial (Eq. 1B), while a second-order multivariable polynomial (Eq. 2B) is represented by the oil recovery objective function. Usually, it is recommended to keep the order of the polynomial model as low as possible, which enhances both the understanding and the prediction of the unknown function.
Therefore, it is important to select the most suitable objective function for achieving the project profitability. While both NPV and oil recovery can be used as objective functions as was previously demonstrated, NPV allows the sensitivity and optimization of parameters that are directly correlated with the profitability of the EOR project, which ensures that the latter is kept within budget. Consequently, NPV as an objective function was further utilized for the sensitivity of the other two sector models (heterogeneous with gravity underride and homogeneous) as well as the subsequent optimization studies.

Heterogeneous model with gravity underride
The same 18 design parameters were considered for the gravity underride model using NPV as the objective function. The RSM was used to build the proxy model for examining the sensitivity of design parameters where about 100 simulation runs were performed. The proxy model generated by RSM algorithm in this case is defined by Eq. (3B) in "Appendix 2," which represents a second-order multivariable polynomial. Oil recovery and related NPV plots were generated, and the proxy model was further verified with most of the points fall closely to the 45-degree line. The tornado chart (Fig. 11) reveals that the three most influential parameters are the tax rate, oil price, and irreducible water saturation. This finding is similar to that obtained for the channeling model. In this case, 13 design parameters out of the 18 parameters are significant. It is worth mentioning that for this model, a design parameter is considered insignificant if its effect on the objective function is within 0.24%.

Homogeneous model
This sector model was chosen as a base model to check the effect of heterogeneity through channeling and gravity underride on the project performance and economics. The same procedure used for the two heterogeneous models was also followed for the homogeneous model. Nevertheless, the heterogeneous design parameters including reservoir heterogeneity (V DP ), cross-flow (K v /K h ratio) as well as transverse and longitudinal dispersivities were excluded. Hence, 14 design parameters were considered for this case.
Oil recovery and corresponding NPV curves were obtained, and the proxy model was verified. This proxy model is given by Eq. (4B) in "Appendix 2." The equation represents a first-order multivariable polynomial that captures the 14 design parameters tested. The sensitivity analysis for the homogeneous model is presented in Fig. 12, which shows that the three most significant parameters are oil price, tax rate, and formation water injected pore volume. In this case, 7 out of the 14 tested design parameters are important. Also, it should be highlighted that for this sector model, a design parameter is considered insignificant if its effect on the objective function is than 1%.

Sensitivity analysis comparison
By comparing the results of the sensitivity analysis study for the three sector models, one can note that there is a consensus on oil price and tax rate being among the most important design parameters (Table 6). Hence, the economic parameters are greatly influencing the objective function (NPV) as opposed to other reservoir and operational parameters. This finding indicates the importance of selecting a proper objective function that truly reflects the profitability of a project. Table 6 shows a comparison between the three sector models with highlighting the top three most influential parameters as well as the parameters influence and percentage effect on the NPV objection function. These most effective parameters include oil price, tax rate, irreducible water saturation, and injected pore volumes of formation water. Oil price is the most influential parameter for the three sector  models with around 60% effect on NPV. This parameter positively affects the NPV where the increase in oil price increases the project revenues and hence increases the NPV. Tax rate negatively affects the NPV where the increase in tax rate increases the costs and hence decreases the NPV. For irreducible water saturation, the decrease in this parameter increases the initial oil saturation, which further increases the size of the prize (amount of hydrocarbons) and hence increases revenues and NPV. For the injected pore volumes of formation water, as this parameter increases, the oil recovery increases and the project revenues increase too, which leads to overall increase in NPV. Nevertheless, field operators should carefully monitor the increase in oil recovery because when it reaches a plateau, the increase in injected PVs of formation water will result in additional costs for water treatment and hence decreases the NPV. Moreover, the results highlight the effect of heterogeneity on NPV through comparing two heterogeneous models (channeling and gravity underride) to a homogeneous model (base model). This is supported by the homogeneous model has the lowest number of influential parameters (7 parameters) compared to the heterogeneous models. The heterogeneous model with gravity underride is the most sensitive model out of the three sector models. This is because this model has the highest number of significant design parameters (13 parameters). Therefore, developing this model requires more attention as opposed to other investigated models, namely channeling and homogeneous models.

Optimization analysis (OPA)
The main point of the experimental design and Response Surface Methodology is to build up the response surface of the objective function which truly represents the response of the objective function. The process variables are divided into state parameters that are usually uncertain and cannot be controlled, and decision variables that are controllable. This step is important as the optimization process depends on the empirical model describing the response variable. Locating the optimal point on the response surface of the objective function is referred as the optimization process (P'azman 1986; Sangvaree 2008; Ghomian 2008; Ghorbani 2008). Schiozer et al. (2019) proposed an approach where the uncontrollable parameters were fixed. The latter approach is more convenient for real case reservoir models as the focus will be on controllable parameters. In this study, the three proposed sector models are synthetic models and hence both uncontrollable and controllable parameters were considered in the process. In this case, the optimization might be used as screening criteria where the results of this study can be applied to any reservoir with similar assumed properties and at certain market conditions. Optimization studies are usually used to produce optimal field development plans and operating conditions, which will either maximize or minimize the objective function selected. In this work, the Designed Exploration and Controlled Evolution (DECE) algorithm was used for the optimization study while considering NPV as an objective function. It is worth mentioning that the optimization study is based on the sensitivity study where only the significant parameters for each model were further utilized.

Heterogeneous model with channeling
As was previously discussed in the Sensitivity analysis (SA) section, 10 sensitive design parameters were considered for optimization of this model (Fig. 9). Figures 13 and 14 show the NPV and oil recovery results of the numerical optimization, respectively. It should be highlighted that 160 simulation runs were performed for this case with the optimal solution in red and the base case in black. Also, the results of the numerical optimization are presented in Table 7.
The findings show a maximum NPV of 1.39 million dollars with a corresponding oil recovery of 79.17% (Figs. 13 and 14). This optimal performance with respect to profitability is best achieved when the EW is injected in the secondary mode, i.e., when the FW pore volumes is zero (minimum) and EW pore volumes is 6.0 (maximum). More so, the corresponding reservoir heterogeneity is 0.6 (minimum), FW residual oil saturation is 0.35, EW residual oil saturation is 0.106 (close to minimum of 0.1), EW sulfate concentration is 21,800 ppm (close to a maximum of 25,000 ppm), irreducible water saturation is 0.35, tax rate is 10% (minimum), oil price is $110 per barrel (maximum), and royalty is 7.94% (close to the minimum of 7.5%). The optimized values for the design parameters using the DECE method physically make sense and are reliable. It is worth mentioning that the insignificant parameters were assigned the base-case values in all optimization studies, as given in Tables 2 and 4 for reservoir and economic inputs, respectively.
Comparing the optimal case with the base case for this channeling model, Fig. 14 shows that the optimal model has higher ultimate oil recovery and faster oil production rate. This difference is basically due to injecting engineered water in a secondary mode as opposed to tertiary injection in the base case. The secondary injection results in an early engineered water response as compared with the tertiary injection where the formation water needs to be swept first before engineered water starts taking effect. This is added to wettability alteration by engineered water, which results in a both favorable fractional flow curve and mobility ratio, and hence better volumetric sweep efficiency at earlier stage. Also, the optimal case can be compared with the base case in terms of economic indicators. Figure 15 shows that the optimal case results in a higher NPV and a shorter payback 1 3 period compared to the base case. The payback period for the optimized scenario is about 2 years as compared to the base case of about 8 years. However, the optimal scenario has a higher maximum exposure compared to the base case, which is due to the high capital and operational costs of applying EWI in the secondary mode during that period.

Heterogeneous model with gravity underride
For the gravity underride model, the screening study results in 13 influential design parameters as previously shown in Fig. 11. Hence, these parameters were further used in the numerical optimization study. The optimization results out of 200 simulation runs are presented in Table 8.
The findings show an optimal NPV of 0.96 million dollars with a corresponding oil recovery of 80. 44%. The optimization shows that for optimal profitability, it is best to inject EW in the secondary, i.e., FW pore volumes is 0.00 (minimum) and EW pore volumes is 5.8 (close to maximum of 6.0). Also, the corresponding EW residual oil saturation is 0.1 (minimum), FW residual oil saturation is 0.3 (minimum), EW sulfate concentration is 23,900 ppm (close to a maximum of 25,000 ppm), reservoir heterogeneity is 0.8, irreducible water saturation is 0.36, selectivity coefficient is 1.0 (maximum), oil price is $110 per barrel (maximum), tax rate is 10% (minimum), royalty is 12.87% (close to minimum), anion exchange capacity is 50, and longitudinal dispersivity is 2.47 (close to minimum).
As expected, the optimal model is more profitable compared to the base. Figure 16 shows that the optimal model has a shorter payback period of 2 years as opposed to 8 years for the base case. The longer payback period for the base case is related to the low oil production rate due to the poor volumetric sweep efficiency where water underrides in the high-permeability layers at the bottom of the reservoir leaving the upper layers poorly swept. Also, the high maximum exposure for the optimal case is due to the large investment on engineered water injected at the start of the project, which increased the costs during this period.

Homogeneous model
For this case, the screening study was performed on 14 design parameters where the heterogeneous parameters were ignored as was previously discussed (Fig. 12). Half of those design parameters were considered for the optimization study where 300 simulation runs were performed. The results of this optimization study are presented in Table 9.
The optimum design has a NPV of 1.44 million dollars with a corresponding oil recovery of 80. 8%. The findings show that the optimal design is achievable when the EW is injected in the secondary mode, i.e., FW pore volumes is zero (minimum). Other optimal design parameters include EW residual oil saturation of 0.1 (minimum), FW residual oil saturation of 0.35 (close to minimum), irreducible water saturation of 0.35, tax rate of 10% (minimum), oil price of 110 dollars per barrel (maximum), and royalty of 8.75% (close to minimum).

3
The optimal design outperforms the base case by achieving faster oil production rate, higher cumulative oil recovery, and shorter project lifetime. Figure 17 supports the previous findings where the optimal design has a higher NPV and shorter payback period of 2.5 years as opposed to that of the base case with 8 years. In terms of maximum exposure, this model is consistent with the previous models where the optimal model has higher maximum exposure compared to the base model due to the costs of early EWI application as was previously discussed.

Optimization analysis comparison
The numerical optimization results are in match with our previous investigations which showed that for best profitability of an EOR project, the EW should be injected in the secondary mode (Adegbite et al. 2017(Adegbite et al. , 2018a. However, higher maximum exposure is expected due to the capital and operational costs related to early EWI application, which can be justified by the corresponding high NPV reward. The latter finding was in match across the three sector models. Moreover, considering the homogeneous model as a base case and comparing both heterogeneous models (channeling and gravity underride), the optimization shows that it is recommended to invest in the channeling reservoir as opposed to the one with gravity underride. This is because the gravity underride is the least profitable reservoir. This finding is supported by the economic indicators summarized in Table 10, especially for the NPV where the channeling reservoir has a value of $1.39 million compared to $0.96 million for that with gravity underride. Furthermore, the channeling reservoir has a lower maximum exposure of −$0.23 million than

Summary and conclusions
Numerical optimization of the engineered water injection (EWI) technique was successfully applied in three sector models (heterogeneous with channeling, heterogeneous with gravity underride, and homogenous) through maximizing the net present value (NPV). Sensitivity analysis was performed on 18 design parameters including reservoir, operational, and economic parameters. This study is one of the very few that discusses the economic aspect of EWI while incorporating the complexity of geochemical reactions through using a mechanistic geochemical model as well as the heterogeneity of carbonates through capturing channeling and gravity underride scenarios. The main findings of this work can be summarized as follows: • The study highlighted that NPV is more representative as an objective function compared to oil recovery where the three optimized models have about similar oil recovery, but different NPVs. • The sensitivity analysis showed that oil price, tax rate, and initial oil saturation are the three most influential design parameters with the most significant effect on the NPV for the three models investigated. • The findings showed that developing the gravity underride model requires more attention as being the most sensitive model with 13 influential design parameters. • The optimization study highlighted the need to apply EWI in the secondary mode of injection in order to achieve the best profitability out of the three models; however, a field operator should be aware of the related high maximum exposure. • High maximum exposure is expected due to the capital and operational costs related to early EWI application, which can be justified by the corresponding high NPV reward. The latter finding was in match across the three sector models.
In the future work, further validation of the developed models can be performed through considering more data in and out of the training set range.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.

Appendix 1
• Marhoun M. A. correlation (1998) to calculate oil formation volume factor at bubble-point pressure: R s is gas solubility (SCF/STB), T is reservoir temperature ( o R), g is specific gravity of solution gas, o is specific gravity of stock-tank oil. • Ahmed T. correlation (1989) to calculate oil formation volume factor at reservoir conditions: B ob is formation volume factor at saturation pressure (bbl/ STB) P is reservoir pressure (psi), P b is bubble-point pressure (psi).

Appendix 2
• Proxy models generated for the three sector models: