Copyright © the author(s). This work is licensed under a Creative Commons Attribution 4.0 International License. DOI: 10.14800/IOGR.429 Received September 5, 2018; revised October 30, 2018; accepted November 8, 2018. *Corresponding author: wanghao1202@petrochina.com.cn 1 Sensitivity Calculation by Adjoint Method: An Application to Reservoir Simulation Hao Wang*, Liaohe Oilfield Company, Petrochina, Panjin, China Abstract Scope of this paper is the sensitivity calculation by using adjoint method. We have a 9×9 two phase two- dimension quarter five spot model with one injector and one producer. The high uncertainty of the permeability field gives us the data mismatch in the model, and we set the permeability as model parameter (uncertainty). The data distribution we get is bottom-hole pressure at injector and water cut at producer for 10 time steps. To get the optimal permeability field, we first calculate the sensitivity coefficients for permeability by using adjoint method. We already have a forward simulator for this problem that is the fully implicit black-oil simulator. Hence, we need to extract necessary information from the forward simulator, i.e. Jacobian matrix, transmissibility, and accumulation etc. Advantage of the adjoint method is that it enables us to reduce the considerable amount of computation time for calculating the sensitivity matrix compared with gradient simulator method. The forward simulation we need at each minimization step is only one time for calculating the sensitivity. Then we minimize the objective function by Levenberg- Marquardt algorithm. Introduction The quality of the reservoir model, i.e., the degree to which it represents the actual reservoir, directly affects reservoir management. This model helps the manager to analyze the behaviors of the reservoir and also it is helpful in production forecasting and optimization. The reservoir model creates on the basis of valuable data obtained during the reservoir life. In the exploration phase, the reservoir model is constructed using 3D seismic data, geologic knowledge of the surrounding area, and log/core measurements from a few exploration wells. These kinds of data are called “static data In the appraisal phase, drilling additional wells provide new information about the reservoir as a well test. Production data and 4D seismic data are available during the production time of a reservoir. These kinds of information are categorized as “dynamic data”. Conditioning reservoir model to the new information obtained about the reservoir is called history matching process, i.e. history matching is the process by which the geological model properties are modified to fit the production data. The objective of history matching is to build a reservoir model that integrates available data and yields production forecasts that are accurate. Incorrectly identifying structural features, such as fluvial channels, can have very serious consequences such as badly placed wells, by-passed hydrocarbon, and failure to find trapped hydrocarbons. A reservoir model is described by many parameters, and each parameter can generate millions of pieces of data. Some parameters are specified per grid block (e.g., permeability and porosity) and others for the entire model or a particular layer (e.g. relative permeability and capillary pressure). Due to insufficiency of available data about the reservoir, history matching is an ill-posed problem. It means that it is possible to obtain reservoir models that honor observed measurements but provide incorrect predictions. To deal with the ill-posed of history matching, the number of parameters has to be reduced. Also, parameterization preserves important geological features and their connectivity that has a significant effect on fluid flow within the reservoir. mailto:wanghao1202@petrochina.com.cn 2 Several methods have been investigated to reduce the number of unknown parameters. Jacquard and Jain (1965) used simple zonation approach to assign one value to a region of the reservoir. Other researchers (Jahns 1966; Bissel et al. 1994; Chavent and Bissell 1998) modified the Jacquard et al.’s method to adaptive one. Some authors (Grimstad et al. 2003; Sahni and Horne 2005) worked on different techniques for parameterization and history matching at different scales. Another powerful approach that is suitable for history matching is KLT. KLT is mathematically defined as an orthogonal linear transformation that transforms a set of possibly correlated data into a smaller number of uncorrelated variables called principal components. KLT is theoretically the optimum transform for given data in least square terms. But for the standard KLT model, it is necessary to carry out an eigen- decomposition of the covariance matrix of the random field, which is expensive for large models. In this research study, we use adjoint method to calculate the sensitivity matrix, thus to reduce the overall computation time during simulation. Figure 1 shows the work flow for the history matching. Figure 1—History matching flow chart. Sensitivity Calculation based on Adjoint method The simulator used in ajoint-based sensitivity calculation is based on fully implicit finite-difference method of the two-phase blackoil model. Several literatures refer to the derivation of adjoint system with fully implicit formulation (Wu et al. 1999; Li et al. 2003). In this chapter, I briefly show the formulation of adjoint method, used in the algorithm for sensitivity calculation. First equation derived from the finite- different equation is as follows. [𝛻𝑦𝑛(𝑓𝑛)𝑇]𝜆𝑛 = −[𝛻𝑦𝑛(𝑓𝑛+1)𝑇]𝜆𝑛+1 − 𝛻𝑦𝑛𝛽,…………..………………..……………….…………(1) where 𝜆𝑛 is the vector of adjoint variables at timestep n, and given by 𝜆𝑛 = [𝜆1 𝑛, 𝜆2 𝑛, … , 𝜆2𝑁 𝑛]𝑇,…………………..………………..……..….……………………………..(2) where 𝑁 is number of gridblock. Returning to Eq. 1, [𝛻𝑦𝑛(𝑓𝑛)𝑇] is transpose of Jacobian matrix at 𝑛, which can be extracted from forward simulation. [𝛻𝑦𝑛(𝑓𝑛+1)𝑇] is the gradient matrix (derivative of finite- difference equation at 𝑛 + 1 with respect to the primary variables at 𝑛). 𝛻𝑦𝑛𝛽 is the sensitivity matrix with respect to the primary variable. Using following initial and boundary conditions, Eq. 1 can be solved backwards in time for n=L, L-1, …, 1. The initial condition and boundary condition are fixed. 𝑑𝑦0 = 0.……….………….……………………………………………………………………………(3) The boundary condition is as follows. 3 𝜆𝐿+1 = 0.……….………………………………………………………………………………………(4) The gradient matrix, [𝛻𝑦𝑛(𝑓𝑛+1)𝑇] is the derivative of accumulation term at 𝑛, because all terms in 𝑓𝑛+1 are independent to the primary variables at 𝑛 except for the accumulation at 𝑛 as shown in Eq. 5. 𝜕𝑓𝑛+1 𝜕𝑦𝑛 = 𝜕{(𝑇𝑡𝑎𝑛𝑠×𝛻𝑝)𝑛+1−(𝐴𝑐𝑐𝑢𝑚𝑛+1−𝐴𝑐𝑐𝑢𝑚𝑛)−(𝑆𝑖𝑛𝑘/𝑆𝑜𝑢𝑟𝑐𝑒) 𝑛+1} 𝜕𝑦𝑛 .........…………………………………..(5) Hence, 𝜕(𝐴𝑐𝑐𝑢𝑚𝑜,𝑖,𝑗 𝑛) 𝜕𝑝𝑖,𝑗 𝑛 = ∆𝑥𝑖,𝑗∆𝑦𝑖,𝑗∆𝑧𝑖,𝑗 1 ∆𝑡𝑛−1 𝑆𝑜,𝑖,𝑗 𝑛 ( 𝜕𝜑𝑖,𝑗 𝑛 𝜕𝑝𝑖,𝑗 𝑛 𝐵𝑜,𝑖,𝑗 𝑛−𝜑𝑖,𝑗 𝑛 𝜕𝐵𝑜,𝑖,𝑗 𝑛 𝜕𝑝𝑖,𝑗 𝑛 𝐵𝑜,𝑖,𝑗 𝑛2 ),…..........….………………………(6) 𝜕(𝐴𝑐𝑐𝑢𝑚𝑤,𝑖,𝑗 𝑛) 𝜕𝑝𝑖,𝑗 𝑛 = ∆𝑥𝑖,𝑗∆𝑦𝑖,𝑗∆𝑧𝑖,𝑗 1 ∆𝑡𝑛−1 𝑆𝑤,𝑖,𝑗 𝑛 ( 𝜕𝜑𝑖,𝑗 𝑛 𝜕𝑝𝑖,𝑗 𝑛 𝐵𝑤,𝑖,𝑗 𝑛−𝜑𝑖,𝑗 𝑛 𝜕𝐵𝑤,𝑖,𝑗 𝑛 𝜕𝑝𝑖,𝑗 𝑛 𝐵𝑤,𝑖,𝑗 𝑛2 ),……….......…….……………….(7) 𝜕(𝐴𝑐𝑐𝑢𝑚𝑜,𝑖,𝑗 𝑛) 𝜕𝑆𝑤,𝑖,𝑗 𝑛 = −∆𝑥𝑖,𝑗∆𝑦𝑖,𝑗∆𝑧𝑖,𝑗 1 ∆𝑡𝑛−1 ( 𝜑𝑖,𝑗 𝑛 𝐵𝑜,𝑖,𝑗 𝑛),…………………......……….………..……………….(8) 𝜕(𝐴𝑐𝑐𝑢𝑚𝑤,𝑖,𝑗 𝑛) 𝜕𝑆𝑤,𝑖,𝑗 𝑛 = ∆𝑥𝑖,𝑗∆𝑦𝑖,𝑗∆𝑧𝑖,𝑗 1 ∆𝑡𝑛−1 ( 𝜑𝑖,𝑗 𝑛 𝐵𝑤,𝑖,𝑗 𝑛).…........………………….………….…………………..(9) The gradient matrix, [𝛻𝑦𝑛(𝑓𝑛+1)𝑇] form the following matrix 𝛻𝑦𝑛(𝑓𝑛+1)𝑇 = 𝜕(𝐴𝑐𝑐𝑢𝑚𝑛)𝑇 𝜕𝑦𝑛 = [ 𝜕(𝐴𝑐𝑐𝑢𝑚𝑜 𝑛)𝑇 𝜕𝑝𝑛 𝜕(𝐴𝑐𝑐𝑢𝑚𝑤 𝑛)𝑇 𝜕𝑝𝑛 𝜕(𝐴𝑐𝑐𝑢𝑚𝑜 𝑛)𝑇 𝜕𝑆𝑤 𝑛 𝜕(𝐴𝑐𝑐𝑢𝑚𝑤 𝑛)𝑇 𝜕𝑆𝑤 𝑛 ].………......….…………………………..(10) Because 𝛽 is assumed as bottomhole pressure at injector and water cut at producer in this project, 𝛻𝑦𝑛𝛽 can be analytically calculated by using following formulation 𝑊𝐶𝑇 = 𝑘𝑟𝑤𝜇𝑜𝐵𝑜 𝑘𝑟𝑜𝜇𝑤𝐵𝑤+𝑘𝑟𝑤𝜇𝑜𝐵𝑜 ,………………………………...…………………………………..………(11) 𝑝𝑤𝑓,𝑖,𝑗 = 𝑝𝑖,𝑗 + 𝑙𝑛(𝑟𝑜,𝑖,𝑗/𝑟𝑤,𝑖,𝑗)+𝑠𝑖,𝑗 (2𝜋)1.127×10−3ℎ𝑘𝑖,𝑗 ( 𝐵𝑚,𝑖,𝑗𝜇𝑚,𝑖,𝑗 𝑘𝑟𝑚,𝑖,𝑗 ) 𝑞𝑚,𝑖,𝑗,……...........………………..…….……………..(12) 𝑟𝑜,𝑖,𝑗 = 0.14036√∆𝑥2 + ∆𝑦2 ..…………………………..…………………………….…………….(13) The adjoint system is ended up computing as shown in Figure 2. The sensitivity coefficients for 𝐽 are given by, 𝛻𝑚𝐽 = 𝛻𝑚𝛽 + ∑ [𝛻𝑚(𝑓𝑛)𝑇](𝜆𝑛)𝐿 𝑛=1 .………………..…………………………………………………(14) We already know 𝜆𝑛 from the computation of Eq. 1. The gradient matrix, [𝛻𝑚(𝑓𝑛)𝑇] is 𝑀 × 2𝑁 sparse matrix as follows 𝛻𝑘(𝑓 𝑛)𝑇 = 𝛻𝑘[(𝑇𝑡𝑎𝑛𝑠 × 𝛻𝑝)𝑛]𝑇 = [ 𝜕𝑓𝑜,1 𝑛 𝜕𝑘1 𝜕𝑓𝑤,1 𝑛 𝜕𝑘1 𝜕𝑓𝑜,1 𝑛 𝜕𝑘2 𝜕𝑓𝑤,1 𝑛 𝜕𝑘2 𝜕𝑓𝑜,2 𝑛 𝜕𝑘1 𝜕𝑓𝑤,2 𝑛 𝜕𝑘1 𝜕𝑓𝑜,2 𝑛 𝜕𝑘2 𝜕𝑓𝑤,2 𝑛 𝜕𝑘2 ⋯ 𝜕𝑓𝑜,𝑁 𝑛 𝜕𝑘1 𝜕𝑓𝑤,𝑁 𝑛 𝜕𝑘1 𝜕𝑓𝑜,𝑁 𝑛 𝜕𝑘2 𝜕𝑓𝑤,𝑁 𝑛 𝜕𝑘2 ⋮ ⋮ ⋮ 𝜕𝑓𝑜,1 𝑛 𝜕𝑘𝑀 𝜕𝑓𝑤,1 𝑛 𝜕𝑘𝑀 𝜕𝑓𝑜,2 𝑛 𝜕𝑘𝑀 𝜕𝑓𝑤,2 𝑛 𝜕𝑘𝑀 ⋯ 𝜕𝑓𝑜,𝑁 𝑛 𝜕𝑘𝑀 𝜕𝑓𝑤,𝑁 𝑛 𝜕𝑘𝑀 ] ........….(15) Figure 3 shows the computation of sensitivity matrix. 4 Figure 2—Computation of adjoint system. Figure 3—Computation of Sensitivity Matrix. Minimization The Levenberg–Marquardt algorithm, which was independently developed by Kenneth Levenberg and Donald Marquardt, provides a numerical solution to the problem of minimizing a nonlinear function. It is fast and has stable convergence. In the artificial neural-networks field, this algorithm is suitable for training small- and medium-sized problems. The Levenberg–Marquardt algorithm blends the steepest descent method and the Gauss–Newton algorithm. Fortunately, it inherits the speed advantage of the Gauss–Newton algorithm and the stability of the steepest descent method. It’s more robust than the Gauss–Newton algorithm, because in many cases it can converge well even if the error surface is much more complex than the quadratic situation. Although the Levenberg Marquardt algorithm tends to be a bit slower than Gauss–Newton algorithm (in convergent situation), it converges much faster than the steepest descent method. The basic idea of the Levenberg– Marquardt algorithm is that it performs a combined training process: around the area with complex curvature, the Levenberg–Marquardt algorithm switches to the steepest descent algorithm, until the local curvature is proper to make a quadratic approximation; then it approximately becomes the Gauss–Newton algorithm, which can speed up the convergence significantly. 5 Gauss–Newton Algorithm. The relationship between Hessian matrix H and Jacobian matrix J can be rewritten as, 𝐻 = 𝐽𝑇𝐽.………………………………………………………………..………………………….….(16) Levenberg's Contribution. Replace the above equation by a "damped version", 𝐻 = 𝐽𝑇𝐽 + 𝜆𝐼,…………………………………………………………………………………………(17) Where λ is always positive, called combination coefficient. I is the identity matrix Levenberg's algorithm has the disadvantage that if the value of damping factor, λ, is large, inverting JTJ + λI is not used at all. Marquardt Modification. Replaced the identity matrix, I, with the diagonal matrix consisting of the diagonal elements of JTJ, resulting in the Levenberg–Marquardt algorithm 𝐻 = 𝐽𝑇𝐽 + 𝜆𝑑𝑖𝑎𝑔(𝐽𝑇𝐽)...…………………………………..………………………………………….(18) From above equation, one may notice that the elements on the main diagonal of the approximated Hessian matrix will be larger than zero. Therefore, with this approximation, it can be sure that matrix H is always invertible. For our problem, we drive the objective function without prior information 𝑂(𝑚) = 1 2 [𝑔(𝑚) − 𝑑𝑜𝑏𝑠] 𝑇𝐶𝐷 −1[𝑔(𝑚) − 𝑑𝑜𝑏𝑠],…..………………………………………….………(19) 𝑂(𝑚 + 𝛿𝑚) = 𝑂(𝑚) + 𝐽𝑇𝑒𝛿𝑚 + 1 2 𝛿𝑚𝑇𝐽𝑇𝐽𝛿𝑚...…………………………………………………….(20) Set 𝜕𝑂(𝑚+𝛿𝑚) 𝜕𝛿𝑚 ≈ 0.…………………………………………………………………………………………(21) We have 𝐽𝑇𝐽𝛿𝑚 = −𝐽𝑇𝑒.………………………………………………………………………………………..(22) Hence, apply Levenberg-Marquardt Algorithm, we get 𝛿𝑚 = −[𝐽𝑇𝐽 + 𝜆𝑑𝑖𝑎𝑔(𝐽𝑇𝐽)] −1(𝐽𝑇𝑒),…………………………………………………………………(23) Where, 𝑂(𝑚) = 1 2 𝑒𝑇𝑒,.......................................................................................................................................(24) 𝑒 = 𝐶𝐷 1/2[𝑔(𝑚) − 𝑑𝑜𝑏𝑠],........................................................................................................................(25) 𝐽 = 𝜕𝑒 𝜕𝑚 ,.................................................................................................................................................(26) 𝐽 = 𝐶𝐷 −1/2 𝜕𝑔(𝑚) 𝜕𝑚 = 𝐶𝐷 −1/2 𝐺𝑙 ,..................................................................................................................(27) 𝐽𝑇𝐽 = 𝐺𝑙 𝑇𝐶𝐷 −1𝐺𝑙 .....................................................................................................................................(28) The simplest way to obtain the correction 𝛿𝑚 is to use Cholesky decomposition on the linear system. The main advantage of the nodal equations is speed. 𝑚𝑙+1 = 𝑚𝑙 + 𝛼𝛿𝑚𝑙+1 ,………………………………………………………………………………..(29) Where, 𝒍 is iteration number. The (non-negative) damping factor, λ, is adjusted at each iteration. 𝑂(𝑚𝑙+1) > 𝑂(𝑚𝑙) → 𝜆𝑙+1 = 𝜌𝜆𝑙,…………………………………………………………………….(30) 𝑂(𝑚𝑙+1) < 𝑂(𝑚𝑙) → 𝜆𝑙+1 = 𝜆𝑙 𝜌 ,………………………………………………………………………(31) 6 where 𝝆 > 1; 𝝀𝟏 is between √ 𝑂(𝑚0) 𝑁𝑑 and 𝑂(𝑚0) 𝑁𝑑 ; f 𝝀 = 0, we have Gauss-Newton Method. If the reduction of objective function is rapid, a smaller value can be used, bringing the algorithm closer to the Gauss-Newton Algorithm, whereas if one iteration gives insufficient reduction in the residual, λ can be increased, giving a step closer to the gradient descent direction. Result Figure 4 and 5 show the WCT and BHP sensitivity comparison between our adjoint sensitivity calculation and sensitivity calculated using perturbation. We are able to capture the trend of both WCT and BHP, the relative difference is generally less than 10%. Figure 4—Sensitivity comparison of WCT. Perturbation Sensitivity 1 2 3 4 5 6 7 8 9 1 2 3 4 5 6 7 8 9 -10 -8 -6 -4 -2 0 2 x 10 -8 Adjoint Sensitivity 1 2 3 4 5 6 7 8 9 1 2 3 4 5 6 7 8 9 -10 -8 -6 -4 -2 0 2 x 10 -8 Sensitivity Difference 1 2 3 4 5 6 7 8 9 1 2 3 4 5 6 7 8 9 -10 -9 -8 -7 -6 -5 -4 -3 -2 -1 0 1 x 10 -9 % Difference 1 2 3 4 5 6 7 8 9 1 2 3 4 5 6 7 8 9 0 0.05 0.1 0.15 0.2 0.25 Perturbation Sensitivity 1 2 3 4 5 6 7 8 9 1 2 3 4 5 6 7 8 9 -10 -8 -6 -4 -2 0 2 x 10 -8 Adjoint Sensitivity 1 2 3 4 5 6 7 8 9 1 2 3 4 5 6 7 8 9 -10 -8 -6 -4 -2 0 2 x 10 -8 Sensitivity Difference 1 2 3 4 5 6 7 8 9 1 2 3 4 5 6 7 8 9 -10 -9 -8 -7 -6 -5 -4 -3 -2 -1 0 1 x 10 -9 % Difference 1 2 3 4 5 6 7 8 9 1 2 3 4 5 6 7 8 9 0 0.05 0.1 0.15 0.2 0.25 Perturbation Sensitivity 1 2 3 4 5 6 7 8 9 1 2 3 4 5 6 7 8 9 -10 -8 -6 -4 -2 0 2 x 10 -8 Adjoint Sensitivity 1 2 3 4 5 6 7 8 9 1 2 3 4 5 6 7 8 9 -10 -8 -6 -4 -2 0 2 x 10 -8 Sensitivity Difference 1 2 3 4 5 6 7 8 9 1 2 3 4 5 6 7 8 9 -10 -9 -8 -7 -6 -5 -4 -3 -2 -1 0 1 x 10 -9 % Difference 1 2 3 4 5 6 7 8 9 1 2 3 4 5 6 7 8 9 0 0.05 0.1 0.15 0.2 0.25 Perturbation Sensitivity 1 2 3 4 5 6 7 8 9 1 2 3 4 5 6 7 8 9 -10 -8 -6 -4 -2 0 2 x 10 -8 Adjoint Sensitivity 1 2 3 4 5 6 7 8 9 1 2 3 4 5 6 7 8 9 -10 -8 -6 -4 -2 0 2 x 10 -8 Sensitivity Difference 1 2 3 4 5 6 7 8 9 1 2 3 4 5 6 7 8 9 -10 -9 -8 -7 -6 -5 -4 -3 -2 -1 0 1 x 10 -9 % Difference 1 2 3 4 5 6 7 8 9 1 2 3 4 5 6 7 8 9 0 0.05 0.1 0.15 0.2 0.25 7 Figure 5—Sensitivity comparison of BHP. Figure 6 and 7 shows the model calibration history and results. We are able to match both BHP and WCT quite well. The updated permeability field keeps the high permeability region in the center. It does not give quite similar results mainly due to the non-uniqueness of solution. Figure 6—History of model calibration (BHP and WCT). Perturbation Sensitivity 1 2 3 4 5 6 7 8 9 1 2 3 4 5 6 7 8 9 -2 -1 0 1 x 10 -4 Adjoint Sensitivity 1 2 3 4 5 6 7 8 9 1 2 3 4 5 6 7 8 9 -2 -1 0 1 x 10 -4 Sensitivity Difference 1 2 3 4 5 6 7 8 9 1 2 3 4 5 6 7 8 9 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 x 10 -5 % Difference 1 2 3 4 5 6 7 8 9 1 2 3 4 5 6 7 8 9 0 0.05 0.1 0.15 0.2 0.25 Perturbation Sensitivity 1 2 3 4 5 6 7 8 9 1 2 3 4 5 6 7 8 9 -2 -1 0 1 x 10 -4 Adjoint Sensitivity 1 2 3 4 5 6 7 8 9 1 2 3 4 5 6 7 8 9 -2 -1 0 1 x 10 -4 Sensitivity Difference 1 2 3 4 5 6 7 8 9 1 2 3 4 5 6 7 8 9 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 x 10 -5 % Difference 1 2 3 4 5 6 7 8 9 1 2 3 4 5 6 7 8 9 0 0.05 0.1 0.15 0.2 0.25 Perturbation Sensitivity 1 2 3 4 5 6 7 8 9 1 2 3 4 5 6 7 8 9 -2 -1 0 1 x 10 -4 Adjoint Sensitivity 1 2 3 4 5 6 7 8 9 1 2 3 4 5 6 7 8 9 -2 -1 0 1 x 10 -4 Sensitivity Difference 1 2 3 4 5 6 7 8 9 1 2 3 4 5 6 7 8 9 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 x 10 -5 % Difference 1 2 3 4 5 6 7 8 9 1 2 3 4 5 6 7 8 9 0 0.05 0.1 0.15 0.2 0.25 Perturbation Sensitivity 1 2 3 4 5 6 7 8 9 1 2 3 4 5 6 7 8 9 -2 -1 0 1 x 10 -4 Adjoint Sensitivity 1 2 3 4 5 6 7 8 9 1 2 3 4 5 6 7 8 9 -2 -1 0 1 x 10 -4 Sensitivity Difference 1 2 3 4 5 6 7 8 9 1 2 3 4 5 6 7 8 9 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 x 10 -5 % Difference 1 2 3 4 5 6 7 8 9 1 2 3 4 5 6 7 8 9 0 0.05 0.1 0.15 0.2 0.25 100 200 300 400 500 600 700 800 900 1000 2975 2980 2985 2990 2995 3000 BHP Time (days) P re s s u re ( p s i) 100 200 300 400 500 600 700 800 900 1000 0 0.2 0.4 0.6 0.8 1 WCT Time (days) Reference Permeability (log10 scale) 1 2 3 4 5 6 7 8 9 1 2 3 4 5 6 7 8 9 Updated Permeability (log10 scale) 1 2 3 4 5 6 7 8 9 1 2 3 4 5 6 7 8 9 0.5 1 1.5 2 2.5 3 3.5 0.5 1 1.5 2 2.5 3 3.5 100 200 300 400 500 600 700 800 900 1000 2975 2980 2985 2990 2995 3000 BHP Time (days) P re s s u re ( p s i) 100 200 300 400 500 600 700 800 900 1000 0 0.2 0.4 0.6 0.8 1 WCT Time (days) Reference Permeability (log10 scale) 1 2 3 4 5 6 7 8 9 1 2 3 4 5 6 7 8 9 Updated Permeability (log10 scale) 1 2 3 4 5 6 7 8 9 1 2 3 4 5 6 7 8 9 0.5 1 1.5 2 2.5 3 3.5 0.5 1 1.5 2 2.5 3 3.5 8 Figure 7—Model calibration results (BHP and WCT). Figure 8 shows the history of objectives. The data misfit decrease dramatically, even in the semilog plot. Figure 8—History of objective (normal and semilog coordinate). Physical Explanation of Sensitivity The effect of permeability change on the BHP and WOR is studied. Figure 9 shows the sensitivity of BHP with the permeability. 100 200 300 400 500 600 700 800 900 1000 2975 2980 2985 2990 2995 3000 BHP Time (days) P re s s u re ( p s i) Reference Initial Updated 100 200 300 400 500 600 700 800 900 1000 0 0.2 0.4 0.6 0.8 1 WCT Time (days) Reference Initial Updated Reference Permeability (log10 scale) 2 4 6 8 1 2 3 4 5 6 7 8 9 0.5 1 1.5 2 2.5 3 3.5 Updated Permeability (log10 scale) 2 4 6 8 1 2 3 4 5 6 7 8 9 0.5 1 1.5 2 2.5 3 3.5 100 200 300 400 500 600 700 800 900 1000 2975 2980 2985 2990 2995 3000 BHP Time (days) P re s s u re ( p s i) Reference Initial Updated 100 200 300 400 500 600 700 800 900 1000 0 0.2 0.4 0.6 0.8 1 WCT Time (days) Reference Initial Updated Reference Permeability (log10 scale) 2 4 6 8 1 2 3 4 5 6 7 8 9 0.5 1 1.5 2 2.5 3 3.5 Updated Permeability (log10 scale) 2 4 6 8 1 2 3 4 5 6 7 8 9 0.5 1 1.5 2 2.5 3 3.5 100 200 300 400 500 600 700 800 900 1000 2975 2980 2985 2990 2995 3000 BHP Time (days) P re s s u re ( p s i) Reference Initial Updated 100 200 300 400 500 600 700 800 900 1000 0 0.2 0.4 0.6 0.8 1 WCT Time (days) Reference Initial Updated Reference Permeability (log10 scale) 2 4 6 8 1 2 3 4 5 6 7 8 9 0.5 1 1.5 2 2.5 3 3.5 Updated Permeability (log10 scale) 2 4 6 8 1 2 3 4 5 6 7 8 9 0.5 1 1.5 2 2.5 3 3.5 0 10 20 30 40 50 60 70 80 90 0 200 400 600 800 1000 1200 1400 Iterations O b je c ti v e OBJ decrease ratio(OBJ initial /OBJ final ) = 0.00099719 0 10 20 30 40 50 60 70 80 90 10 0 10 1 10 2 10 3 10 4 Iterations O b je c ti v e OBJ decrease ratio(OBJ initial /OBJ final ) = 0.00099719 9 Figure 9—Sensitivity of BHP. We can see from the colorbar that most the grid cells have negative value, which means if the permeability increases, the BHP decreases. It is easy to understand that higher permeability is more conductive for water than lower permeability. Besides, with the time increasing, the cells around the producer are becoming more sensitive to the BHP. Figure 10 shows the sensitivity of WOR with the permeability. Figure 10—Sensitivity of WOR. 10 We can see there is a channel connecting injector and producer. The sensitivity of these channel cells are positive, which means if the permeability increases in these areas, the WOR increases. More water tends to flow in this channel. The sensitivity of cells in the edge part is negative, meaning higher permeability in these areas can decrease the WOR. More water flow to this area. Effect of α and λ in the LM minimization process. α is the searching step in the LM algorithm. Small α results in more iterations while large α may lead to convergence failure. Large searching step may make the derivative stepping cross the minimization point. Figure 11 shows the comparison of history matching process between α=0.01 and 3. Figure 11—Comparison of history matching process between α=0.01 and 3. The left part of Figure 11 shows the history matching process with α equals 0.01. We can see the converging process is slow. The blue line is the final curve. It is limited by the iteration numbers. If given sufficient iterations, it can converge. The right part shows history matching process with α equals 3. We can see the matching curve jumps around, and the convergence is failed. Figure 12 below shows the objective function behavior with iteration numbers for nine different α. We can see if α equals 0.2, the objective function converges slowly and smoothly, although 68 iterations are needed. But if the α equals 3, the objective function fluctuates without tendency to converge. These are two extreme cases. Additionally, if α equals 1.7, it converges very fast with only 8 iterations. Figure 12—Objective function behavior with iteration numbers for different α. 11 λ in the LM algorithm is the perturbation added to the original function, to avoid the singular matrix. If λ is too small, the original function may still be singular, while if λ is big, more iterations will be needed to converge the objective function. Figure 13 below shows the trend. Figure 13—Objective function behavior with iterations for different λ. We can see if the λ equals 0.001, the objective function fluctuates a lot, and the code gives the warning that the objective matrix is badly conditioned. The result is not correct. For λ of 10, it converges slowly and smooth, but 72 iterations are needed. For λ of 20, it does not converge with the limited iterations, but it will converge given sufficient iterations. Additionally, λ of 0.5 gives very fast convergence. We can see both α and λ affect the history matching process, and good combination of them is desired. Conclusions 1. We successfully incorporate the sensitivity calculation with adjoint method into the software, and the sensitivity results are quite close to the results calculated by perturbation. 2. The adjoint method only need one-time simulation to give the sensitivity calculation, instead of M times (the number of parameters), which show great advantages in computational efficiency. 3. We can carry out history matching using Levenberg-Marquardt method to match both WCT and BHP. We are also able to keep the high permeability trend of the permeability field. 4. More suitable search step lengths and perturbation in Levenberg-Marquardt method will give faster converge. Integrated with adjoint method, it is a fast and accurate history matching workflow. Conflicts of Interest The author(s) declare that they have no conflicting interests. References Wu, Z., Reynolds, A.C., and Oliver, D.S. 1999. Conditioning Geostatistical Models to Two-Phase Production Data. Paper presented at SPE Annual Technical Conference and Exhibition, New Orleans, Louisiana, September 27-30. SPE-56855-MS. Li, R., Reynolds, A.C., and Oliver, D.S. 2003. History Matching of Three-Phase Flow Production Data. SPE Journal 8(4):25-31. SPE-87336-PA. Oliver, D.S., Reynolds, A.C., and Liu, N. 2008. Inverse Theory for Petroleum Reservoir Characterization and History Matching: Front Matter. Journal of Affective Disorders 2008(175):488-493. Jacquard, P. 1965. Permeability Distribution from Field Pressure Data. SPE Journal 5(4): 281-294. SPE-1307-PA. 12 Jahns, O.H. 1966. A Rapid Method for Obtaining a Two Dimensional Reservoir Description from Well Pressure Response Data. SPE Journal 6(4): 315-327.SPE-1473-PA Chavent, G. and Bissell, R. 1998. Indicator for the Refinement of Parameterization in Inverse Problems in Engineering Mechanics. Ed. M. Tanaka and G.S. Dulikravich, 309-314. Oxford, UK: Elsevier Science. Grimstad, A.A., Mannseth, T., Nævdal, G., et al. 2003. Adaptive Multiscale Permeability Estimation. Computational Geosciences 7(1): 1-25. Sahni, I. and Horne, R.N. 2005. Multiresolution Wavelet Analysis for Improved Reservoir Description. SPE Reservoir Evaluation & Engineering 8(1): 53-69. SPE-87820-PA. Hao Wang is a petroleum engineer at Liaohe Oilfield Company, a subsidiary of PetroChina. Wang’s specialties include viscous oil development technologies, like SAGD and steam drive. Wang holds a bachelor’s degree in petroleum engineering from China University of Petroleum, Beijing, and master’s degree in petroleum engineering from Northeast Petroleum University.