Sensitivity Analysis and Optimization of Surfactant-Polymer Flooding under Uncertainties

— Analyse de sensibilite´ et optimisation sous incertitudes de proce´de´s EOR de type surfactant-polyme`re — La re´cupe´ration ame´liore´e des hydrocarbures, base´e sur l’injection d’agents chimiques, est actuellement une solution conside´re´e comme une des plus prometteuses pour ame´liorer la re´cupe´ration de champs matures. Pour des proce´de´s base´s sur une injection de polyme`res et de surfactants, plusieurs parame`tres doivent eˆtre pris en compte pour pouvoir estimer le retour sur investissement : concentrations des espe`ces chimiques injecte´es, tailles des bouchons, saturation en huile re´siduelle, taux d’adsorption des agents sur la roche, etc. Parmi ces parame`tres, certains sont des parame`tres de controˆle et d’autres des parame`tres incertains. Pour des ope´rateurs, de´terminer les valeurs optimales des parame`tres de controˆle tout en tenant compte des incertitudes lie´es aux parame`tres incertains n’est pas une taˆche facile dans la pratique. Cet article propose une me´thodologie


INTRODUCTION
Enhanced Oil Recovery (EOR), especially chemical EOR, is currently one of the most active sector of research in the oil industry.Different techniques based on chemicals injections such as polymers, surfactants with or without alkalines can in some cases help increase significantly the recovery rate of mature reservoirs, where classical production techniques like water injection are no more effective.In this domain, numerical modeling can be used at different stages.At the beginning of an EOR project, screening models based on simplified assumptions [1] or comprehensive reservoir models [2] can be used to appraise if one of the previous techniques may be appropriate for a given reservoir.Such models can then be employed to probe different injection scenarii and to optimize them in order to maximize the production or the project rentability.
Specifically, the key features of chemical EOR are mainly three-fold: -lower the interfacial tension in order to mobilize the residual oil by means of surfactant, -limit the surfactant consumption/adsorption by means of alkali and/or dedicated surfactants, -improve the mobility ratio/sweep efficiency by means of polymer.Like for other classical production techniques used in the oil industry, the optimization of chemical processes is not easy in practice because of a large number of uncertain parameters.These parameters can be related to rock properties (porosity, absolute or relative permeabilities, etc.), in-situ fluid properties (water salinity, viscosities, densities, etc.), economic conditions (capital expenditure, facility costs, oil price, etc.), operating conditions (maximum injection rate, limiting bottom-hole pressure, etc.).Moreover, when injecting polymers, surfactants or alkalines, the adsorption of these chemical species on the rock, the Capillary Desaturation Curve (CDC) or the variations of the oil-water interfacial tension with respect to salinity and surfactant concentration, for instance, are poorly known at the reservoir scale.
Despite these uncertainties, which usually come from a lack of data, measurements accuracy or scaling from core to field scale, operators have to determine the most appropriate design paramaters (type of chemicals, concentrations, initiation times, slug sizes, etc.) for a given EOR process and a given reservoir.Optimization workflows based on Response Surface Modeling (RSM) can help choose an optimal design.Different methodologies have already been proposed in previous works [3][4][5][6] for EOR processes based on P, SP or ASP injections of polymer and/or surfactants with or without alkalines (A stands for alkaline, S for surfactant and P for polymer).These workflows consist of three main steps.First, an experimental design (Box-Behnken, three-level D-optimal experimental designs, etc.) is defined to sample the values of the uncertain and design parameters.From the simulation results of the experimental design, a quadratic Response Surface (RS) is built to approximate the production results of a reservoir simulator with respect to the input parameters.At last, the RS are used to find the optimal values of the design parameters making the uncertain ones vary.The use of such surfaces enables to reduce the computational cost since reservoir simulations are no more performed at the final step.More complex than quadratic RS, kriging or Gaussian process models, were used in other reservoir contexts [7] or in other applications [8].Reference [7], for instance, uses Gaussian processes to approximate key outputs of numerical models from a fixed number of simulations.Adaptive design strategies have also been proposed instead of fixed designs to improve Gaussian process approximations for given responses [9].
In this work, a RSM technique is used to perform the sensitivity analysis of a synthetic but realistic reservoir model where slugs of polymer and surfactant are injected.The parameters adopted for this analysis are the concentrations of chemicals, their adsorption rates on the rock, the CDC and, more generally, the relative permeability curves.At this stage, a particular parameterization is introduced in order to perform the sensitivity analysis with respect to these properties.Indeed, pointwise modifications would require too many parameters and make the RSM methodology useless.On that occasion, a new way of approximating the production curves by response surfaces is proposed.Instead of computing response surfaces for each time of well production data, a Karhunen-Loe`ve decomposition is applied on the curve and response surfaces are used to approximate the coefficients of this decomposition.In our case, this leads to a significant reduction of the computation time.Let us note that this methodology combining Principal Component Analysis (PCA) and Gaussian process modelling of PCA coefficients has been recently proposed by [10] and applied in order to perform the sensitivity analysis of an atmospheric dispersion model.
In a second part, the response surfaces obtained in the previous step are used to optimize the concentrations to be injected by taking economic constraints and uncertainties into account.
In the following of this paper, we introduce our base model in Section 1.The definition and the construction of the response surfaces are described in Section 2 and the results obtained in our example are also presented in this section.In Section 3, we detail our methodology to perform optimization under uncertainties and introduce an economic model to illustrate how an optimal design based on RSM can be obtained.Concluding remarks and future works are summarized at the end.The distributions of porosities and permeabilities are depicted in Figures 1 and 2. The permeability map is a log-normal realization ranging from 1.7 to 5 710 mD with a mean of 265 mD.The porosities vary from 1.7 to 13.5% with an average value of 7.6%.Figure 1 also displays the locations of the wells which are placed according to a direct five-spot pattern with one injector in the middle of the reservoir and four producers in the corners.A dead-oil model is used to represent the fluid properties.Oil and water viscosities are equal to 2 and 0.43 cP and oil and water densities are equal to 0:824 g=cm 3 and 1 g=cm 3 (in-situ conditions).
Initially, the water saturation is equal to the residual water saturation S wi ¼ 0:2 everywhere.Oil is produced by means of three injection sequences over a 2 000-days production period: -a water flooding up to day 506 (% 0.4 pore volume (PV)); -a slug of surfactant and polymer is injected from day 506 to day 886 (% 0.3 PV); -a second slug containing only polymer is injected from day 886 to day 2 000 (% 0.86 PV).
A constant injection rate of 150 m 3 =d is imposed at reservoir conditions and, for each producer, a maximum production rate of 37:5 m 3 =d is fixed at reservoir conditions along with a minimum bottom-hole pressure of 15 bar.
We decide to stop the water flooding at day 506 because, at that time, the water-cut at well P 4 reaches 90% (see Fig. 9), whereas the field water-cut is 50%.
All simulations were performed using the reservoir simulator PumaFlow TM .For the case considered in this work, one single simulation runs in a few minutes.But the simulation times become much longer when working with real data defined on larger grids with several hundreds of thousands of grid blocks.For these cases, the simulation times can take a few hours.
In the following, we briefly describe how the surfactant-polymer injection is modeled within PumaFlow TM in order to introduce our uncertain and design parameters afterwards.For more details about this model, the reader can refer to [11,12].
In chemical EOR, surfactants are usually injected along with water to reduce the interfacial tension between oil and water.Polymers can also be injected at the same time or afterwards.They intend to increase the water viscosity and thus to reduce the mobility ratio to improve the sweep efficiency of the water front.Polymer and surfactant also contribute to reduce the residual oil saturation S or after water flooding.
In our example, the mobility reduction R m due to the polymer is defined as a function depending on the polymer concentration C pol and the oil-water interfacial tension r ow as a function of the surfactant concentration C surf .Both functions are represented in Figures 3 and 4.
The reduction of residual oil saturation S or is modelled by means of a Capillary Desaturation Curve (CDC) describing the decrease of the residual oil saturation with respect to the capillary number N c , defined as the ratio of the viscous driving force over the interfacial tension: In this definition, l w is the water viscosity without polymer, R m is the reduction of the water mobility induced by the polymer, jũ w j jjis the norm of the water phase velocity and r ow is the oil-water interfacial tension.In our specific case, we assume that the polymer only modifies the water visocity and does not alter the relative permeabilities; therefore R m only increases the water viscosity.
We assume that the CDC can be represented using the following relation: where erfcðxÞ ¼ 2 ffiffi p p R þ1 x e Ày 2 dy is the complementary error function, S or;w the residual oil saturation after water flooding without polymer and surfactant and S or;c the residual oil saturation after chemical flooding.Note that other CDC parameterizations have already been proposed and the reader can refer to [13] for other examples.The CDC in our initial model is represented in Figure 5, where S or;w ¼ 0:32, S or;c ¼ 0, a ¼ 0:875 and N c0 ¼ 5 Â 10 À3 .Here, we assume that the irreducible water saturation S wi remains constant during surfactant-polymer flooding.
Another consequence of the reduction of the oilwater interfacial tension by surfactants is the modification of the relative permeability curves.At very low interfacial tension values, oil and water phases should be miscible and their relative permeabilities tend to be cross-shaped.To mimic this effect, we model the relative permeabilities by means of Corey power laws and make vary the Corey exponents for two different capillary number values N c1 and N c2 .Specifically, for N c < N c1 , a Corey exponent n C;w ¼ 2 is used for the water phase and a Corey exponent In our example, as illustrated in Figure 5, N c1 and N c2 are automatically computed in a preprocessing step before the reservoir simulation, according to the locations of the two levels S or;w and S or;c in the CDC curve.For N c1 N c N c2 , the reservoir simulator deduces the shape of the oil and water relative permeability curves in each gridblock by interpolating the curves of these two sets according to the local value of N c .
We also assume that the maximum water relative permeability kr w;max depends on N c taking: where erf ðxÞ ¼ 1 À erfcðxÞ is the error function, kr w;max;w the maximum water relative permeability during a simple water flooding and kr w;max;c the maximum water relative permeability during chemical flooding.
Polymers and surfactants may be partly adsorbed on the rock, which tends to reduce their efficiency in practice.In our case, the polymer adsorption rate A pol;max is taken constant and equal to 50 lg=g.The surfactant adsorption rate A surf is modelled as a function of the surfactant concentration using a Langmuir isotherm given by: with A surf ;max the maximum adsorption rate, A surf ;L the Langmuir coefficient, and C surf ;m the mass concentration of surfactant.In the initial model, we choose A surf ;max ¼ 500 g and A surf ;L ¼ 0:1 m 3 =kg.Because of a high permeability area linking the injector and the fourth producer, the water breakthrough first occurs at this well.The effect of the first slug of surfactant and polymer noticeably appears for the wells P 3 and P 4 where it manages to reduce the rise of the water rate.But the injected concentrations are clearly not sufficient here.Figure 10 shows the recovery obtained with the initial model with and without the injection of chemicals.The overall volume of oil in place being equal to 154 126 m 3 , the recovery factor obtained after day 2 000 is 56.2% using only water and rises by 0.8% with chemicals.

First Numerical Results
In the next section, we make the concentrations vary along with some potentially uncertain parameters in order to observe the variability of the production.A response surface is then built in order to find out the most influential parameters.

SENSITIVITY ANALYSIS
Starting from our base model, we now assume that some parameters are uncertain or need to be optimized.These parameters are: -the maximum polymer and surfactant adsorption rates A pol;max and A surf ;max (see Eq. 3); -the coefficient a which drives S or and kr w;max curvatures (see Eq. 1, 2);  -the capillary number N c0 where S or and kr w;max reach their mean values; -the water and oil phases Corey exponents n C;w;c and n C;o;c for N c ! N c2 ; -the injected surfactant and polymer concentrations C surf and C pol;sp during the second injection sequence; -the injected polymer concentration C pol;p in the third and last slug.Table 1 gives the range of variations of these parameters along with their value in the base model.
In practice, A pol;max , A surf ;max , a, N c0 , n C;w;c and n C;o;c are difficult to quantify whereas C surf , C pol;sp and C pol;p are design parameters to optimize in order to obtain the most profitable production scheme.When working on real field cases, this list of parameters is far from being exhaustive.For example, slug sizes and initiation times are design parameters which are also often taken into account to optimize a process [6].Here, we restrict our list of uncertain and design parameters to Table 1.In fact, the instrinsic limitation of Gaussian process regression with standard kernels lies in the number of parameters it can handle.Recent developments on kernels based on an analysis of variance (ANOVA) decomposition and component selection have been proposed [14] and may be a way to extend our list of parameters.

Variance-Based Sensitivity Analysis and Gaussian Process Approximation
Sensitivity analysis aims at identifying which inputs most influence the output of a given model.
S i (also called primary effect of X i ) is the first order Sobol index which measures the part of the output variance due to X i only.For i 6 ¼ j, S ij measures the part of the output variance due to the interaction between X i and X j .The indices of higher orders S ijk are defined in the same way.Sobol indices are included in the interval ½0; 1 and their sum is equal to one when all input variables are independent.A parameter X i or the interaction between two parameters X i and X j are all the more influential on the output variance that their Sobol indice is close to one.The total sensitivity index S Ti can also be used to measure the overall sensitivity of the response with respect to a parameter X i .This index is defined by: where k#i denotes all terms where index i appears.Different techniques (FAST, quasi Monte-Carlo, etc. [16]) exist to estimate the first and total indices.However, all these estimation methods usually require several thousands of calls to the model: when it is expensive, the Sobol indices can not be estimated directly.In practice, a standard way to perform this task is to build an approximation of the model output with a limited number of model evaluations, and to use this cheap approximation in all sensitivity calculations.In the literature, such an approximation is called a proxy, a responsesurface model, an emulator, or a meta-model.Among all proxy models, a popular method is the Gaussian process regression, also referred to as kriging [17].Let us denote by yðxÞ any scalar output of the reservoir simulator (e.g. for a given well, the cumulative oil at a given time).In this setting, yðxÞ is considered to be a realization of the Gaussian process: where: h stands for the mean of Y which is often modelled as: b j h j ðxÞ h j being given functions; -ZðxÞ is a stationary Gaussian process with zero mean and covariance function covðZðx 1 Þ; Zðx 2 ÞÞ.In practice, the covariance function is specified with a parametric form.For example the power-exponential covariance function is given by: In the previous equation, x with: For more details, the reader can refer to [17].Estimates of b and r 2 are given by: At the end, the value of the response surface for the input parameters x Ã is given by RSðx Ã Þ ¼ EðY ðx Ã ÞjSAÞ and the conditional variance error can be used as the error of this predictor.
Practically, the experimental design is obtained with an optimal Latin Hypercube technique with maximum property [18,19].It is well known that this type of space-filling design is suited for exploration but not for hyperparameter estimation of the Gaussian process model.Recent attempts aim at combining space-filling designs and designs dedicated to hyperparameter estimation.But no satisfying solution exists so far and this point is beyond the scope of this paper.
The overall predictability of the response surface is usually measured by the so-called Q 2 coefficient.If an independent test sample ðx j ; y j Þ j¼1;...;n test is available, the Q 2 coefficient can be estimated by: where RS denotes the response surface built using the current training sample.Alternatively, another estimation method consists in using a cross-validation technique (cv): where RS Àj is the response surface built from the current training sample deprived of the point ðx j ; y j Þ.If the predictability is too low with the initial training sample, it is possible to adopt a sequential strategy: additional points are simulated in regions where the predictor is not accurate until the overall Q 2 is large enough.Generally, a Q 2 over 0:9 ensures a quite good quality of the response surface.
Coming back to our problem, we want to estimate the sensitivity indices of the recovery at several times with respect to some parameters.This means that we need to build one different proxy for each time t i with i ¼ 1; . . .; m, where an output yðt i ; xÞ is computed by the reservoir simulator.However, the computation time related to these calculations increases with the size of the sampling n (see Eq. 4) and also with m since an optimization of the hyper-parameters is performed at each time.To reduce drastically the computation times related to m, one solution is to characterize the time response yðt; xÞ by its most relevant components using a Karhunen-Loe`ve (KL) decomposition [20].In this type of decomposition, the output yðt; xÞ is expanded in the following way: where the / k basic functions are orthonormal and sorted in descending order according to the explained variance of the output.In practice, this decomposition is computed with a standard PCA algorithm.
In our approach, we apply the RS not to each yðt i ; xÞ but to each component a k ðxÞ of the decomposition.In practice, only a few components are needed for an accurate reconstruction, meaning that we only have to build a small number of proxy models (e.g.k ¼ 1; :::; 10).The results presented in what follows were obtained using this approach.

Application
In this section, our RSM methodology is applied to the oil recovery of our five-spot case.A classical LHS composed of 400 simulations was performed.Figure 11 shows the recovery factors computed during this experimental design.
A KL decomposition where only the first two components were kept, was then carried out.In our case, these two components are sufficient to explain more than 99% of the variance of the reponses of the initial LHS.This number of components should be compared to the number of time steps which here amounts to 202.A speed up of 101 is thus obtained and, because of the size of our LHS, the construction time of the RS based on the KL decomposition turns out to be much faster than with the classical approach.For both approaches, attempts were made to use smaller LHS with 90 and 100 simulations.Unfortunately, the quality of the responses was not satisfactory, which compelled us to increase the LHS significantly.Note that an adaptive experimental design was not used here but that it could help reduce that dimension.For more details on that topic, the reader can refer to [9].
To validate our response surfaces, 192 confirmation runs were used on the whole: two LHS of 90 and 100 simulations, the case of a pure water flooding and another case known to give a good recovery.The Q 2 test coefficients obtained with these runs are given in Figure 12.During the first water flooding up to day 506, the response does not vary and Q 2 test is set to 1 by convention.When the injection of the first slug starts, this value decreases suddenly because of the low value Response variance Time (days) Q 2 test coefficients (in blue) and variance (in red) of the responses obtained from 192 confirmation runs.
of the variance but it is not problematic because the simulator response does not change a lot in that part (see definition of Q 2 test given by Eq. 6).From day 1 000, all Q 2 test coefficients are greater than 0.9.Figures 13 and 14 show the productions estimated by the response surface and the reservoir simulator for two confirmation runs.Note that Q 2 test values equal to 0:9157 and 0:9090 were obtained for the two components a 1 and a 2 of the KL decomposition.Associated cross plots of the 192 confirmation runs are given in Figures 15 and 16.
The primary and total effects of each parameter were computed from day 1 000 using RBD-FAST technique [21].Their evolution is shown in Figures 17 and 18. Globally, the three concentrations, the maximum polymer adsorption and the capillary number N c0 have the most important primary effects.When considering the total effects, we notice that this hierarchy is still the same.But values which sometimes exceed the primary effects by 0:1 suggest that interactions between parameters also have an influence on the responses.The calculations of the secondary effects revealed that the parameters with the highest primary effects interact with each other and that A surf ;max , a, n C;w;c , n C;o;c are not influent at all in our case.Table 2 shows the highest interactions of order 2. Practically, the computation of the Sobol indices can give hints to reduce the number of uncertain parameters before performing an history-matching or an optimization (see Sect. 3).As an illustration, the shape of the response surface at the final time is given in Figures 19 and 20.In the first one, each plot gives the variation of the surface with respect to one parameter setting the other ones to their mean values.The second figure shows the variations of the surface with respect to both A pol;max and C surf keeping the other parameters constant.

OPTIMIZATION UNDER UNCERTAINTIES
In this last part, we illustrate how to optimize our SP process using the RSM achieved in the previous section.
Starting from the total production obtained with a pure water flooding, we estimate thanks to the RS the additional production we can obtain with our SP process.This difference of production is then used by a simplified economic model to compute the Net Present Value (NPV).The parameters of this model are given in Table 3.Here, two scenarios with an oil price equal to 10 $=bbl and 100 $=bbl are considered.For both scenarios, we compute the distribution of the NPV and the distributions of the optimal concentrations C surf , C pol;sp , C pol;p taking the uncertainties related to A pol;max , A surf ;max , a, N c0 , n C;w;c , n C;o;c into account.
The results are shown in Figures 21-28.They were obtained by taking a random sample of 500 values of the uncertain parameters and by maximizing the NPV for each value of the sampling.Note that for these optimizations: -the uncertain parameters as well as C surf , C pol;sp , C pol;p were constrained to the ranges given in Table 1; -uniform laws were used for the uncertain parameters; -a BFGS algorithm was used to compute the optimal concentrations for each sample of the uncertain parameters.
For an oil price equal to 10 $=bbl, the SP process proposed in our example is obviously not profitable.It becomes much more cost-effective in the second scenario.In that case, the RSM suggests to inject the maximum polymer concentration and a intermediate value around 4 000 ppm for the surfactant.

CONCLUSIONS
In this work, a methodology was proposed to perform the sensitivity analysis and the optimization under uncertainties of various parameters related to SP processes.The list of the parameters considered in this paper is far from being exhaustive and also depends on the simulator chosen for a field study.However we saw, through our example, how suitable parameterizations can allow a reservoir engineer to carry out a sensitivity analysis and an optimization with functional inputs like relative permeabilities.On the other hand, we proposed a methodol-ogy of RSM which consists in a KL decomposition of the simulator time response and in an approximation of the components of this decomposition by a Gaussian Variations of the response surface at the final time with respect to each parameter.
Variations of the response surface at the final time with respect to A pol;max and C surf .
process.Compared to more classical approaches which aim at modeling several different times of the response independently, this approach can reduce the computation times when trying to approximate the whole reservoir output.We also showed that a good predictability of the response surface could be achieved.But our example revealed that the simulator response (here the cumulative oil production) with respect to the chosen uncertain and design parameters can be quite complex and that the number of simulations required to obtain a good predictability was high compared to the number of parameters.As a following of this work, adaptive strategies could be investigated to reduce the size of the experimental design.Finally, we saw how to make use of the computed response surfaces to optimize control parameters like polymer and surfactant concentrations.Beyond the conceptual nature of our test case, it would be interesting to extend the set of parameters (e.g. chemical flooding initiation time, salinity-dependent surfactant adsorption and/or interfacial tension) and to apply this workflow to real data and to other numerical models.
To illustrate our methodology, we consider, throughout this work, a synthetic two-dimensional reservoir model.Its dimensions are 505 m Â 505 m Â 10 m and it is discretized into 101 Â 101 Â 1 gridblocks of constant sizes equal to 5 m Â 5 m Â 10 m.The reservoir top depth is constant with z top ¼ 2 095 m.

Figure 1
Figure 1 Porosity distribution and well locations.

Figure 4
Figure 3Mobility reduction as a function of polymer concentration.

Figures 6 -Figure 5 Figure 6
show the oil-and water-cut obtained with our initial model at the four producing wells.The starts of the second and third injection sequences, where first a slug of surfactant and polymer and then a slug of polymer are injected, are also displayed in the figures.

Figure 7 InitialFigure 8 Figure 9
Figure 7Initial model oil-and water-cut at producer P 2 .
Figure 11Recovery factors computed during the experimental design.
Figure 13Comparison of the 25 th confirmation run of the first LHS (90 points) with the response surface.
Figure 15Simulated (a Sim 1 ) and predicted values (a GP 1 ) of a 1 for 192 confirmation runs.

Figure 17 Evolution
Figure 18Evolution of the total Sobol indices.

FigureFigure
FigureDistribution of C surf (ppm) with an oil price equal to 100 $=bbl.
sensitivity indices, which evaluate the impact of the inputs on the variance of the model output.More precisely, let us assume that the inputs X are modelled as random variables with a given probability distribution and denote Y the model output which, by propagation, is also a random variable.The variancebased sensitivity indices, or Sobol indices [15], are defined by: Several measures have been proposed in the literature, ranging from derivative-based indicators to norms of conditional probability density functions.Here, we focus on variance-based