Screening method using the derivative-based global sensitivity indices with application to reservoir simulator

Résumé — Méthode de criblage basée sur les indices de sensibilité DGSM : Application au simulateur de réservoir — Les simulateurs de réservoir peuvent impliquer un grand nombre de paramètres d’entrées. Cependant, il est possible de réduire cette complexité en ne se focalisant que sur les paramètres inﬂuents. La détection de ces paramètres peut être e ﬀ ectuée grâce à l’analyse de sensibilité. Néanmoins, pour des problèmes de grande dimension, les méthodes classiques d’estimation des indices de sensibilité nécessitent un très grand nombre d’évaluation du simulateur de réservoir, ce qui est très couteux en temps de calcul. Récemment, pour pallier au problème de dimensionnalité, de nouveaux indices de sensibilité ont été déﬁnis. Ces indices dit DGSM (pour derivative-based global sensitivity measures) sont basés sur la moyenne des dérivées partielles et ont un lien avec les indices totaux qui sont des indices de sensibilité basés sur la décomposition de la variance. Cet article introduit une version révisée des indices DGSM, ce qui permet de proposer une méthode de criblage e ﬃ cace et moins couteuse en temps de calcul que les méthodes classiques. Aﬁn d’apprécier l’e ﬃ cacité de la méthode proposée, des résultats sur des fonctions analytiques ainsi que sur un cas synthétique de réservoir sont présentées. Abstract — Screening method using the derivative-based global sensitivity indices with application to reservoir simulator — Reservoir simulator can involve a large number of uncertain input parameters. Sensitivity analysis can help reservoir engineers focusing on

Abstract -Screening method using the derivative-based global sensitivity indices with application to reservoir simulator -Reservoir simulator can involve a large number of uncertain input parameters. Sensitivity analysis can help reservoir engineers focusing on the inputs whose uncertainties have an impact on the model output, which allows reducing the complexity of the model. There are several ways to define the sensitivity indices. One of the quantitative definition is the variance based sensitivity indices which can quantify the amount of output uncertainty due to the uncertainty of inputs. However, the classical methods to estimate such sensitivity indices in a high dimensional problem can require a huge number of reservoir model evaluation. Recently, new sensitivity indices based on averaging local derivatives of the model output over the inputs domain have been introduced. These so-called derivative-based global sensitivity measures (DGSM) have been proposed to overcome the problem of dimensionality and are linked to total effect indices which are variance based sensitivity indices. In this work we propose a screening method based on a revised DGSM indices which increase the interpretability in some complex case and have a lower computational cost, as demonstrated by numerical test cases and by an application on a synthetic reservoir test model.

INTRODUCTION
Reservoir simulators are complex computer codes that model the physical laws governing the recovery process, and which are mainly modelled by mathematical equations for the three phases flow (oil, gas and water) through porous media. These simulators involve a large number of input parameters. The information gathered on such inputs comes from direct measurements, which are clearly very limited and are marred by considerable uncertainty. Thus, it is important to detect influential inputs, whose uncertainties have an impact on the model output. Once identified, one can reduce the complexity of the model by fixing the non-influential inputs on default values (defined by experts) and focus the attention on the influential inputs. Sensitivity analysis (SA) is the study of how the variation (uncertainty) in the output of the computer model can be apportioned, qualitatively or quantitatively, to different sources of variation in the input of the model. Put in another way, it is a technique for systematically changing parameters in a model to determine the effects of such changes on the output. The local SA methods refer to the study of the sensitivity at a fixed point in the input domain, typically the simple derivative ∂Y/∂x (i) of the output Y with respect to a given input X (i) taken at some fixed point x 0 in the input domain. The global SA (GSA) methods [17][18][19]26] refers to the sampling-based methods in which the model is evaluated for combinations of values sampled from the distribution (assumed known) of the inputs. Once the sample is generated, several strategies (including simple input-output scatterplots) can be used to derive global sensitivity measures for the factors. The variance-based methods are nonlinear with respect to the input parameters, and are based on analysis of variance (ANOVA) decomposition which is the decomposition of the total variance V of output into terms due to individual factors plus terms due to interaction among inputs. Most variance-based methods are quantitative, and in this work we will focus on this class of methods, and more specifically on Sobol's indices. One of the main issues with variance based methods is computational time. Indeed, a reservoir simulator is often very costly in terms of computational time. Furthermore such variance based methods require generally several thousands simulations that are usually not affordable in common applications. In order to perform SA with a limited number of runs, metamodel methods can then be used. In the latter the simulator input/output relation is approximated using different statistical regression techniques starting from an initial set of carefully chosen training runs. Then, if a reasonably good approximation is obtained, the estimated metamodel is used instead of the complex simulator to compute the sensitivity indices. Metamodel methods have known a quick development in the last decade and different approaches have been suggested in many different scientific disciplines [2,3,10,16,20,[27][28][29][30]. However, despite significant advances in the area, construction of a sufficiently accurate approximation for high dimensional computer code using relatively low number of model evaluation is problematic.
Screening methods aim at reducing the input dimensionality by identifying the non-influential inputs with a low computational cost in terms of model evaluation. Screening design proposed by Morris [11] is adapted for high dimensional expensive computer model. This method is an one factor at time (OAT) technique that vary one input parameter at a time and measures the impact on the output. Indeed, the method is based on calculating a sensitivity index called elementary effect, which provide a good compromise between accuracy and efficiency. However, even if this method is computationally cheaper than other SA methods, it involves hundreds or thousands (depending on the number of inputs and the complexity of the model) of model evaluation, which is still computationally intensive with realistic reservoir simulators for which each simulation requires several hours or days.
Recently, Sobol and Kucherenko [24 25] have proposed new sensitivity indices based on averaging local derivatives of the model output over the inputs domain. It was shown that the so-called Derivative-based Global Sensitivity Measures (DGSM) can be easily estimated and much faster than the global sensitivity indices. In addition, different methods exist to efficiently compute the derivatives of reservoir simulators, in this case DGSM represent a valid alternative to the Morris method for screening the input parameters. In this work we propose a revised derivative-based sensitivity index that allows a better convergence of the estimation and increases the interpretability in some complex cases. We propose a screening method based on the defined indices. We then employ the method to perform a screening of a high-dimensional numerical test cases and a reservoir simulator.

GLOBAL SENSITIVITY ANALYSIS AND MORRIS DESIGN
First let consider a mathematical model for a reservoir simulator where Y is a scalar output of the computer code, X = (X (1) , . . . , X (d) ) a unit d-dimensional input vector (X ∈ [0, 1] d ) which represents the uncertain parameters/factors of the simulator and f : [0, 1] d → R is a function that models the relationship between the input factors and the output of the computer code.

Global sensitivity analysis
The main idea from Sobol approach [26] is to decompose the response Y = f (X) into summands of different dimensions via analysis of variance decomposition (ANOVA) defined as where f 0 is a constant, f j 's are univariate functions representing the main effects, f jl 's are bivariate functions representing the two way interactions, and so on. The integrals of every summand of this decomposition over any of its own variables is assumed to be equal to zero, i.e.
Using the orthogonality, Sobol [26] showed that such decomposition of f (X) is unique and that all the terms in (2) can be evaluated via multidimensional integrals where E(Y) and E(Y|X ( j) ) are respectively the expectation and the conditional expectation of the output Y. Analogous formulae can be obtained for the higher-order terms. If all the input factors are mutually independent, the ANOVA decomposition is valid for any distribution function of the X (i) 's and using this fact, squaring and integrating (2) over [0, 1] d , and by (4), we obtain where is the variance of the conditional expectation that measures the main effect of X j on Y and Variance-based sensitivity indices, also called Sobol indices, are then defined by where 1 ≤ j 1 < . . . < j s ≤ d and s = 1, . . . , d. Thus, S j = V j /V is called the first order sensitivity index (or the main effect) for factor X ( j) , which measures the main effect of X ( j) on the output Y, the second order index S jl = V jl /V, for j l, is called the second order sensitivity index expresses the sensitivity of the model to the interaction between variables X (i) and X ( j) on Y and so on for higher orders effects. The decomposition in (8) has the useful property that all sensitivity indices sum up to one.
The total sensitivity index (or total effect) of a given factor is defined as the sum of all the sensitivity indices involving the factor in question.
where #i represents all the S j 1 ,..., j s terms that include the index i. Total effect index of an input X (i) measures the part of output variance explained by all the effects in which it plays a role. Note however that the sum of all S T i is higher than one because interaction terms are counted several times. It is also important to note that total effect indices can be computed by a single multidimensional integration and do not require computing all high order indices (see Sobol [26]). Then comparing the total effect indices provides information about influential parameters. Indeed, one can suppose that the input is non-influential if its total effect S T i is less than 0.01. GSA enables to explain the variability of the output response as a function of the input parameters through the definition of total and partial sensitivity indices. The computation of these indices involves the computation of several multidimensional integrals that are estimated by Monte Carlo method and thus requires huge random samples. For this reason GSA techniques are prohibitive if used directly using the computer code (fluid flow simulator for example).

Morris's screening method
The screening method introduced by Morris [11] is based on a one factor at time (OAT) experimental design. The points of the Morris design are sampled from a d-dimensional plevel grid, as the range of each input X (i) is divided into p equal level. The impact of varying one input at a time is evaluated by the so-called elementary effect that, for i = 1, . . . , d, is defined as where ∆ is a multiple of 1/(p − 1) with p the number of levels, X r is a randomly chosen point in [0, 1] d such that The group of points composed of X r and X r i 's are called trajectories. Thus, the Morris design is structured in R random trajectories composed of R(d + 1) points. The sensitivity measures proposed by Morris Morris [11] are defined as a statistics of the elementary effect. The first one is the meanμ iμ which is a measure of the ith input importance. The second statistic is the standard deviation of the elementary effectσ î which is a measure of the non-linearity and the involved interactions of the ith input. However,σ i does not allow to distinguish between non-linearities and interactions.
Noting that when the model is non-monotonic the elementary effects of opposite signs cancel each other, Campolongo et al. [4] proposed the sensitivity measureμ * i , which is a revised version ofμ iμ * To identify the non-influential inputs, the sensitivity measuresμ * i andσ i are simultaneously considered.

Derivative-based global sensitivity measures
First introduced by Sobol and Gresham [21] and then studied in Kucherenko et al. [9], Sobol and Kucherenko [24 25], DGSM are a new sensitivity indices based on averaging local derivatives of the model output over the inputs domain.
Assume that ∂ f (X)/∂x (i) , for i = 1, · · · , d, are squaredifferentiable. The DGSM indices are defined as Thus, calculation of DGSM indices is based on the evaluation of integrals, which is easily performed using a classical Monte Carlo (MC), quasi-Monte Carlo (QMC) or latin hyper-cube (LHS) sampling. The empirical estimator of ν i is given byν

Link between DGSM and GSA
Recently, Sobol and Kucherenko [24] have established the link between the DGSM index ν i and the total effect index S T i for input variables following uniform and normal distributions. Here, we assume that . . , d, the link between V T i and ν i is defined by the following inequality Thereby the total effect indices have the following upper bound where V is the total variance of the model. If Υ i 0, then X (i) can be considered as non-influent input, which make Υ i a good candidate for a screening procedure.

Link between DGSM and Morris method
Kucherenko et al. [9] have introduced two derivative-based sensitivity indices that are very similar to the Morris indices µ i and σ i . These indices are defined, for i = 1, . . . , d, as whereM i is equivalent to µ i andΣ 2 i equivalent to σ 2 i . In addition, these indices are more accurate than Morris's indices, which cannot correctly consider effects with characteristic dimensions less than ∆. Indeed, in (13) the elementary effects d r i are calculated as finite differences with the increment δ, which has the same order of magnitude as the uncertainty range of inputs. In contrast with the derivativebased indicesM i andΣ 2 i , where the elementary effects are substitute by the local derivatives. We can also note that the DGSM indices ν i can be defined as

Refining the DGSM index Υ
In the previous section we have shown that the sensitivity index Υ is an upper bound of the total effect index. However, for some complex model Υ can be much larger than the corresponding total effect index. In this case, it is difficult to decide which inputs are influential and which are not. In addition, Υ estimation involves the variance of the model output V. Empirical results (see next section) show that estimation of V requires more model evaluation than estimation Initialization -Build an initial design X using QMC -Run the reservoir simulator at X and store the output data and the gradients calculation We propose here a different version of Υ that we call Υ * , which is defined, for i = 1, . . . , d, as This index is a normalized upper bound of V T i . Indeed, the link between Υ * i and V T i is defined as In addition, Υ * has the following useful properties The drawback of Υ * is the loss of the link with the total effect indices. Nevertheless, the use of Υ * offers a stronger measure to define the non-influential inputs, as we will see in the next sections.

DGSM-BASED SCREENING METHOD
As shown in section 3, estimating Υ * for each input allows one to detect the influential inputs in the model. Indeed, one can state that ifΥ * i ≤ 0.01, the corresponding input X (i) can be defined as non-influential one. However, such criteria may be too strong in the case of very high dimensional model (typically more than 100 input parameters), the fact is that the sum of effects due to inputs with small sensitivity indices may be significant on the output model variance. Because of that, and by using the property (27) of the Υ * measures, it is more robust to state that the influential inputs are the set D of d * inputs whose the Υ * i are the most high and which respect the following criteria i∈D Υ * i ≈ 0.98 (28) Given that the reservoir simulator evaluation may be computationally demanding, it is important to use a sequential strategy to build the design of point by reusing at each step the already evaluated points. The use of the so-called QMC Sobol sequence [22,23] is an efficient way to build a sequential design. Our choice is motivated by two main properties of the QMC Sobol sequences. First, this technique is based on the generation of deterministic quasi-random sequences with a good space-filling property of the unit hypercube, in other words the input domain are well covered for fairly small sets. Second, the points of the Sobol sequence are independent. That is, by enriching the design sequentially, one keeps the space filling properties of the Sobol sequence. Since a sequential method computes successive estimation of the Υ * i indices, a practical test is needed to determine when to stop the iteration. In this work we propose to use the following error criteria where vectors Υ * (l) = (Υ * 1 , . . . , Υ * d ) are the lth estimation of Υ * i indices and · is the Euclidean norm. Thus, we define the stopping criteria as err ≤ 0.05. A schematic representation of the entire screening method is shown in Figure 1. Note that one can use the same stopping criteria (29) to estimate Υ i indices.

NUMERICAL TESTS
In this section, two numerical test cases are used to demonstrate the estimation performance of the DGSM index Υ *

The test case of Morris
The test function proposed by Morris [11] contains 20 input parameters and defined as follow β i = 20 f or i = 1, · · · , 10 β i, j = −15 f or i, j = 1, · · · , 6 β i, j,l = −10 f or i, j, l = 1, · · · , 5 β i, j,l,s = 5 f or i, j, l, s = 1, · · · , 4 the remaining β i and β i, j are independently generated from a standard normal distribution. The remaining β i, j,l and β i, j,l,s are sets to zero. In Figure 2 one can see the results of computing Υ i and Υ * i sequentially with a QMC design ranging from n = 5 to n = 256. In addition, the sample size when the stopping criteria of the proposed screening method is reached is represented by the red vertical line. Thereby, we can notice that Υ * i converge faster than Υ i . Indeed, the stopping criteria is reached at n = 43 for Υ * i and at n = 67 for Υ i estimations. In Table 1, for the inputs selected by the proposed screening method, the values of the total effect indices (obtained by the so-called extended-FAST method [18] using a sample of size N = 3.5 × 10 4 ) as well as the values of Υ i and Υ * i (computed at the stopping criteria sample size n = 67 and n = 43) are reported. Notice that the values of the indices S T 11 , . . . , S T 20 are smaller than 0.005 and then the corresponding inputs are considered to be non-influential. It can be seen that for this test case both indices (Υ i and Υ * i ) are able to identify the influential inputs correctly at the stopping criteria (29). Furthermore, even if at n = 67 the indices Υ i are under estimated they are almost for all inputs greater but close to the total effect indices. To conclude on this numerical test, we can say that for the Morris function the developed screening method based on Υ * i is efficient and not consuming in terms of model evaluations. However, estimating Υ i indices provides more information since it is very close as upper bound to the total effect indices but the drawback is the additional model evaluation cost.

The g-Sobol function
Consider the g-Sobol function, which is strongly nonlinear and is described by a non-monotonic relationship. Because of its complexity and the availability of analytical sensitivity indices, this function is a well-known test case in the studies of GSA. Let define the g-Sobol function for 200 input parameters as follow to the variability of the model output is represented by the weighting coefficient a k . The lower this coefficient a k , the more significant the variable X (k) .
The analytical values of Sobol's indices are given by where 1 ≤ j 1 < . . . < j s ≤ d and s = 1, . . . , d. The analytical values of the total effect indices are shown in table (2). Figure 3 shows the sequential estimation of Υ * i and Υ i indices with a QMC design ranging from n = 5 to n = 256. As for the previous test example the estimation of Υ * i indices converge faster than those Υ i . In this test case the stopping criteria (29) is reached at n = 41 for Υ * i and n = 59 for Υ i . The analytical values of the total effect indices as well as the values of Υ i and Υ * i (estimated at the sample size n = 59 and n = 41) for the inputs which are identified as influential by the screening method are reported on Table 2. Notice that the values of the indices S T 12 , . . . , S T 200 are smaller than 0.007 and then the corresponding inputs are considered to be non-influential. One can see that the informations provided by the indices Υ 1 , · · · , Υ 4 are difficult to interpret. Indeed, for the g-Sobol function the values of Υ i 's provide only a qualitative informations, because for some inputs Υ i > 1 which is higher than the maximal value for the total effect indices. These results may due to the model non linearity with respect to the inputs. On the other side, despite the non-linearity and non-monotonicity of the model the Υ * i measures performs very well in terms of quantitative interpretability. The above two numerical tests, show us that both DGSM indices are adapted to identify the non-influential inputs. Moreover, the (27) property of Υ * i indices allows one to use an automatic screening method regardless the complexity of the studied model.

RESERVOIR FORECASTING APPLICATION
In this section, the proposed screening method is applied on a reservoir simulator. As the goal here is to apply the method on a high dimensional case, we choose to use horizontal and vertical permeability as inputs parameters. However, since the number of grid blocks in the considered reservoir simulation model is large, we applied the most basic parametrization technique which is the zonation to reduce the dimension of the problem. This technique consists in dividing the reservoir into relatively small number of zones (subregions)  and assume that each zone is homogeneous. In other words, one fixes the permeability (horizontal or vertical) over all the grid blocks of the considered zone.

Reservoir model description
The PUNQS case is a synthetic reservoir model taken from a real field located in the North Sea. The PUNQS test case, which is qualified as a small-size model, is frequently used as a benchmark reservoir engineering model for uncertainty analysis and for history-matching studies Floris et al. [7]. The geological model contains 19 × 28 × 5 grid blocks, 1761 of which are active. The reservoir is surrounded by a strong aquifer in the North and the West, and is bounded to the East and South by a fault (see Figure 4). A small gas cap is located in the centre of the dome shaped structure. The geological model consists of five independent layers, where the porosity distribution in each layer was modelled by geostatistical simulation. The layers 1, 3, 4 and 5 are assumed to be of good quality, while the layer 2 is of poorer quality. The field contains six production wells located around the gas-oil contact. Due to the strong aquifer, no injection wells are required. For more detailed description on the PUNQS model, see PUNQS [14].
As input parameters, we considered the horizontal and vertical permeability of 60 zones (12 for each layer). Thus, we have a model of 180 inputs, which are supposed independent and defined as follow -Z1, · · · , Z60: horizontal permeability in X direction -Z61, · · · , Z120: horizontal permeability in Y direction -Z121, · · · , Z180: vertical permeability The values of the permeability at each zone are distributed uniformly over where PZ i is the arithmetic mean of the permeability values of the grid blocks which compose the ith zone. The analyzed output is the production watercut data (the proportion of water in the produced oil) after 20 years of production of the well 5 which the perforation location correspond to the grid blocks [17; 11; 3 : 4], where the notation 3 : 4 means that there is a perforation at layer 3 and at layer 4 . The reservoir test model was run using the PumaFlow T M [13] simulator, which allows one to compute gradients using gradient simulator method [1,8] with an additional ≈ 33% of the simulation time per each calculated gradient.

Screening
For this reservoir model, we first performed a convergence study of the DGSM indices. Υ i and Υ * i are computed sequentially with a QMC design ranging from n = 5 to n = 200. From Figure 5 one can see that Υ * i indices converge much faster than Υ i indices. The stopping criteria (29) was reached at n = 18 for Υ * i and at n = 42 for Υ i indices estimation. The screening method identifies 13 zones as influential. The DGSM sensitivity measures of these 14 parameters are reported in Table 3. We can see that these zones correspond to the region where the studied well is located and the closest north region. This result demonstrates the relevance of application on a reservoir simulator of the developed screening methodology. To corroborate these screening results we build a metamodel using a standard implementation of Gaussian process method (GP). The GP code used here is a commercial version implemented in the CougarFLow T M software [5]. For more detail on the technical aspect of the used GP we refer to section 3 of Busby et al. [3]. The GP metamodel is build using the results obtained on the QMC design of the size n = 200. However, instead of using the full design of 180 inputs we just selected the 13 inputs identified as influential by the screening method. Thus, rather than building metamodel that approximate a reservoir model of 180 inputs we build a metamodel f which involves only the parameters that supposed to be influential on the output. To assess the prediction accuracy of the metamodel we performed an extra 100 random evaluations of PUNQS simulator (with 180 inputs) and compare the simulator results to the metamodel ones. The measure of the accuracy is given by the Q2 criteria defined by where y i denotes the ith simulator evaluation on test set,ȳ is their empirical mean and f (x i ) is the predicted value at the design point ). The empirical Q2 criteria of the considered metamodel f is equal to 0.94, which means that the metamodel explain 94% of the output variance. Thus the obtained metamodel is sufficiently accurate to perform a global sensitivity analysis. In the second column of Table 3 the reported total effect indices was computed through the metamodel f and using extended-FAST method. From Table 3 one can say that the developed screening method permits to detect the most important inputs. Concerning the inputs Z106 whose total effect index (estimated using the metamodel) is equal to zero, we can suspect that the influence of this parameter has been underestimated by the considered metamodel. We can notice that for this reservoir model the values of Υ i and Υ * i indices have the same magnitude for all selected inputs at the stopping criteria and the values of Υ i are smaller than the estimated total effect which is due to the underestimation of Υ i indices at n = 42.

Input
Total effect   Figure 5: Convergence of the Υ i and Υ * i indices estimates versus the sample size for the production watercut output after 20 years of production of the well 5. Graphics of the third line correspond to the Convergence of the Υ i and Υ * i indices estimates of the non-influential inputs

CONCLUSIONS AND DISCUSSIONS
In this work, we presented a new sequential screening method which is based on DGSM indices. We defined a new DGSM index Υ * i in order to have a stronger quantitative measure to define the non-influential inputs. We also used the QMC Sobol sequence sampling method, which allows an intelligent sequential estimation of the DGSM indices in order to reduce the number of model evaluation. We empirically showed, by applying on two analytical models and a reservoir synthetic test case, that the proposed screening method is efficient to detect the non-influential inputs for an acceptable computational cost. Computing DGSM indices requires model gradients estimation. A classical way to compute the derivatives is to use the finite-difference approximation method. However, this method suffer from the fact that the required number of model evaluation is equal to n(d + 1), where d is the number of inputs and n the number of points where derivatives are estimated. Since a reservoir simulator evaluation is gener-ally time consuming, the finite-difference method is infeasible for models with a high number of inputs (roughly more than 20). Therefore, it is clear that the ability to calculate derivatives efficiently is important to estimating DGSM indices within an acceptable computational cost. In the framework of reservoir simulation different methods have been developed for more or less computationally efficient gradients calculation. In this paper we utilized a reservoir simulator, which allows to compute gradients using the direct method, or also called in reservoir engineering the gradient simulator [1,8]. This method is based on the solution of the governing analytical finite difference equations of flow, which automatically calculates the gradients during the simulation with an additional ≈ 33% of the simulation time per each calculated gradient. Thus the required number of model evaluation to estimate DGSM indices is ≈ n(d + 1)/3. However, the most efficient method to calculate the gradient of a functional with respect to the reservoir simulator parameters is the adjoint state method when this functional depends on those models parameters through state variables, which are the solution of the differential equation that define the problem. The advantage of this method comparing to the gradient simulator method is that it consists of the computation of one unique extra linear system and the computation of the gradient with respect to the model parameters is equivalent to one evaluation of the simulator. In other words, it means that the computational cost of gradient calculation is independent of the number of model parameters. So the required number of model evaluation to estimate DGSM indices is equal to 2n. For more details on the mathematical aspect of the adjoint state method and its applications in reservoir simulation we refer to [6,12,15]. In addition to further testing on reservoir models, using a reservoir simulator which allows to compute gradients using adjoint state method is a topic of futur work.