|Fluid Dynamics & Materials Processing|
An Improved Parameter Dimensionality Reduction Approach Based on a Fast Marching Method for Automatic History Matching
1Southern Marine Science and Engineering Guangdong Laboratory, Zhanjiang Bay Laboratory, Zhanjiang, 524000, China
2Zhanjiang Branch of China National Offshore Oil Corporation, Zhanjiang, 524000, China
3College of Petroleum Engineering, Yangtze University, Wuhan, 430100, China
*Corresponding Author: Deng Liu. Email: firstname.lastname@example.org
Received: 25 September 2021; Accepted: 26 October 2021
Abstract: History matching is a critical step in reservoir numerical simulation algorithms. It is typically hindered by difficulties associated with the high-dimensionality of the problem and the gradient calculation approach. Here, a multi-step solving method is proposed by which, first, a Fast marching method (FMM) is used to calculate the pressure propagation time and determine the single-well sensitive area. Second, a mathematical model for history matching is implemented using a Bayesian framework. Third, an effective decomposition strategy is adopted for parameter dimensionality reduction. Finally, a localization matrix is constructed based on the single-well sensitive area data to modify the gradient of the objective function. This method has been verified through a water drive conceptual example and a real field case. The results have shown that the proposed method can generate more accurate gradient information and predictions compared to the traditional analytical gradient methods and other gradient-free algorithms.
Keywords: History matching; parameter dimensionality reduction; sensitive area; gradient correction
Reliable reservoir models can accurately reproduce the complete reservoir developing process and provide convincing guidance for production management decision-making [1–5]. History matching that can correct reservoir model parameters based on geological information and production data is a critical step in reservoir numerical simulation process. Traditional history matching methods are generally time-consuming, and the accuracy could hardly be guaranteed. In recent decades, automatic history matching methods has made great progress with improved optimization algorithms, in which case the model parameters would change automatically to minimize the difference between the established model and the real one. Currently, automatic history matching technology has become a hot research area in intelligent reservoir development.
Generally speaking, there is a close correlation between reservoir physical properties and well production performance, and history matching algorithms are mainly used to calculate the gradient and to direct the countermeasures. To the best of our knowledge, the widely used algorithms for history matching mainly include gradient-class method [6,7], set-class algorithm such as EnKF , ES-MDA  and SPSA , interpolation algorithms such as NEWUOA  and Wedge , and global optimization algorithms such as PSO  and GA . The gradient-class method shows very high convergence rate, in which case, calculate the gradient of objective function or Hessian matrix is the key. However, there are also two limitations when using the gradient-class method for history matching: (1) gradient solving is difficult. The reservoir model usually contains hundreds of thousands of unknown parameters, which makes it impractical to store and calculate the sensitivity matrix. (2) pseudo correlation between real gradient and calculated gradient. The prior reservoir model based on scarce geological information is generally quite different from the real reservoir, which leads to inaccurate covariance between calculated production data and reservoir parameters. That means the calculated correlation between distant grids is large instead of small. As a result, the calculated gradient deviates seriously from the actual gradient.
In order to reduce the dimensionality of variables in the history matching, several parameterization methods were proposed. The overall idea is to convert the original high dimensional optimization problem to a low dimensional problem without changing the main features. The parameterization methods mainly include principal component analysis (PCA) method [15–18], Karhunen-Loeve (K-L) decomposition method , gradual deformation method (GDM) , pilot point method , discrete cosine transformation (DCT) method [22–24] and so on. PCA is a statistical analysis method, which is widely used in areas such as finance , machine learning  and so on. This method adopts orthogonal transformation to linearly transform the observed values to a series of possibly related variables, so as to project them into the values of a series of linearly uncorrelated variables. These unrelated variables are called principal components. In the oilfield, PCA method is commonly applied to extract reservoir heterogeneity and to reduce the calculation dimension.
The pseudo correlation between gradients can be weakened by dividing the single-well sensitivity area and then modifying the calculated gradient . Single-well sensitive area refers to the grid area that would affect the well production performance, and the grid outside the area seldomly has any impacts. FMM is an effective approach to determine single-well sensitivity areas [28–30]. In FMM, the Eikonal equation is established in accordance with the static reservoir parameters. The propagation time of pressure wave from the well node to different grids is calculated to track the boundary migration. And then the critical area (sensitive area) affects the well production performance is obtained. During the production history matching, only the reservoir parameters in the single-well sensitive area need to be revised.
This paper presented a two-step history matching method by integrating the FMM and the PCA. The brief calculation scheme of history matching is shown in Fig. 1. In this method, FMM is introduced to localize the reservoir to distinguish the single-well sensitivity area. Subsequently, based on multiple prior geological models, PCA is used to parameterize the reservoir model and to reduce the dimensionality. The gradients of the parameterized variables is then modified based on the information on sensitivity area. The objectives of this paper are, first, apply the FMM method to distinguish the single-well sensitive areas. Second, to establish a mathematical model for history matching and to propose a correction method to correct the gradient. Third, to verify the feasibility of the method with a water drive conceptual example and a real field case.
2.1 FMM Method
The pressure wave diffusion equation (Eq. (1)) was derived by using the asymptotic ray theory of geometrical optics and seismology .
where, k is permeability, is fluid viscosity, P is pressure, is porosity, Ct is the comprehensive compressibility coefficient, t is the time variable. Fourier transform and asymptotic expansion of Eq. (1) can prove that under the high frequency limit, and the propagation of pressure wave front can be described by the Eikonal equation (Eq. (2)) [28–31]:
where, T is the diffusion flight time at the x position; is the diffusion factor; is the unit conversion factor with a constant value of 0.0853. FMM tracks the leading edge of pressure wave by solving the Eikonal equation. Herein, upwind difference method is introduced to solve Eq. (2):
where, T is the flight time at position (i, j, k); Δx, Δy and Δz are the grid size on the x, y, z directions, respectively; vx, vy and vz are the propagation velocity in the x, y, z directions, respectively; Tx, Ty and Tz are the grid tracking time in the x, y, z directions, respectively. As shown in Fig. 2, by considering reservoir permeability anisotropy in Cartesian grid system, the process to track the flight time is as following:
(1) Set the initial start grid and then mark it as a freezing grid (red point in Fig. 1a);
(2) Search for adjacent grids and calculate the flight time T;
(3) Find the grid to which the minimum flight time is required (“1” point in Fig. 1b and mark it as new freezing grid);
(4) Search for adjacent grids surround the freezing grid. Then calculate the flight time and set a new freezing grid (“2” point in Fig. 1d);
(5) Repeat Steps 2 to 4 until all grids are marked as freezing grids.
In order to reduce the calculation cost, we set the time threshold of single-well sensitive area to a preliminarily defined tracking range. In this paper, we arranged 1/T it from large to small and calculated the cumulative value. When the cumulative value accounted for 99% of the sum, all the grids with values less than the time threshold are set as the initial tracking range of the well.
2.2 Principle of History Matching
In this section, we first briefly described the mathematical model of history matching and the process to do parameter dimensionality reduction. Then, we explained how to use the single-well sensitive area data to correct the gradient of the objective function.
2.2.1 Mathematical Model of History Matching
History matching is a typical inverse-problem solving process, which aims to generate the maximum posteriori estimate of the reservoir model by matching the observation data. The correction between the actual observation data and the model parameters is:
where, dobs is the observed data, such as water content, oil production rate, bottom hole pressure, etc.; g( ) refers to the numerical simulation process; m is the model parameter; ed is the measuring error. According to the Bayesian framework, the conditional PDF of m conditional dobs can be written as:
where, CM is the correlation matrix of model parameters; CD is the covariance matrix of measuring errors. Thus, in the history matching process, the following formula should be minimized:
where, mprior is the prior estimation vector of reservoir model. In the history matching process, the prior model was often used as the initial value of the continuous iteratives. Since mprior often follows the multivariate Gaussian distribution in practical applications, the average values of these prior model parameter that represented the prior geological features would be optimized as the initial value. In this paper, is used instead of mprior:
In general, the actual reservoir parameters often have high dimensionality, making it pretty difficult to directly calculate the CM. In the past, many methods such as K-L decomposition method, discrete cosine transform method and gradual deformation method have been applied to reduce the dimensionality of the parameters in history matching. In this paper, efficient SVD decomposition method in PCA was introduced to parameterize the reservoir model and to deconstruct the CM to make Nw non-zero singular values. CM was calculated as following:
where, is the diagonal matrix, in which diagonal elements are singular values; is the singular matrix. Herein, we defined a vector w, and m could be expressed by w:
The objective function can be approximately converted to:
The maximum posteriori estimation mMAP could be inversely calculated after the maximum posteriori estimation wMAP was obtained by minimizing Eq. (11):
2.2.2 Correction Method to Correct Gradients
In the numerical simulation process, the production data generated by the reservoir model based on the model parameters m generally followed the following relationship:
where, g(m) is the simulated data; G is the sensitivity matrix of g(m). The gradient of Eq. (11) can be written as
where, Gw is the sensitivity coefficient matrix of g(w). At the l iteration step, the variable could be updated by:
where, αl is the step size in search. In the previous study, the sensitivity matrix G was constructed over the m field. After parameterization, it was necessary to transfer the sensitivity matrix over m field to over w field, that was, transferring G to Gw. Herein, we defined ρ(m) as the sensitivity localization function over the m-field: . The gradient of ρ(m) is:
where, Nd is the number of the observation data. After the single-well sensitive area was obtained by FMM, the sensitivity gradient beyond the sensitivity area was defined as . Based on the flight time, a exponential model was used to obtain the single-well sensitivity localization correlation model inside the sensitivity area:
where, Fj is the sensitive area of well j generated by FMM; Tref is the maximum flight time within the sensitive area of well j, which is related to the static parameter field m. It was also worth noting that the flight time could be normalized to make the localization model more applicable.
Then we defined ρ(w) as the sensitivity localization function over the w field. Set:
According to the chain rule to solve the partial derivatives:
where, . Therefore, localized analytical gradient considering the single-well sensitive area could be obtained by the Schur product between and . The latest updated variables were:
3 Case Study
3.1 Conceptual Case
A two-dimensional heterogeneous reservoir model with three reservoir fluids, oil, gas and water was first constructed to verify the proposed method. The size of this model was set to be 20 × 30 × 1. The grid size is, Dx = Dy = 80 m, and Dz = 30 m. We selected twenty prior reservoir realizations and a real model that provided the observation data. The real permeability field of the model was shown in Fig. 4a, which was of strong heterogeneity with high permeability zones. Fig. 4b was the average permeability field of twenty reservoir models. The average model was to be used as the initial model for history matching. Fig. 4c was the permeability field of a specific model. For all of models, the average porosity was 0.2, the initial water saturation was 0.2, and the gas was dissolved gas. More fluid physical properties could be seen in Table 1 and Fig. 3. There were one injection well and five production wells. The injection rate was a constant of 500 m3/d. The liquid flow rate of the production wells, 100 m3/d. A total of 1200 parameters including permeability and porosity were to be estimated by matching the actual production data.
In the history matching process, the FMM method was firstly used to calculate the flight time. Taking the injection well in the average model as an example (Fig. 5), the flight time was found to be related to permeability field. According to Eq. (1), the calculated migration velocity of the pressure wave was also related to the reservoir permeability. Generally speaking, higher permeability corresponds to better reservoir physical properties, and the flight time would be shorter. Based on the flight time, the sensitive area of each single well was then determined to localize the reservoir, as shown in Fig. 6. It was also worth mentioning that it only took 1.3 s to calculate the flight time without any numerical simulation runs, indicating the higher efficiency of FMM for sensitive area division as compared to the traditional methods based on experience, seismic data  or Voronoi diagram .
Subsequently, the PCA method was used to reduce the dimensionality of the parameters included in the reservoir model. In order to verify the effectiveness of the proposed method to weaken the pseudo-correlation, the finite difference method (FDM), simultaneous perturbation stochastic approximation (SPSA) [34,35], localized analytic gradient method (LAGM) and analytical gradient method without considering sensitive area (AGM) were adopted respectively to estimate the gradient of the objective function and to perform the history matching with 70 iterations. AGM refers to that the variable is updated by Eq. (15), and LAGM refers to that the variable is updated by Eq. (21). Some details about SPSA algorithm were listed in the Appendix A. For the four algorithms, the iteration step was set to 0.5. Figs. 7 and 8 were the calculated gradient by the four algorithms and its gradient correlation coefficient of each variable at the 10th and the 70th iteration step. Among them, the gradient calculated by FDM marked by red circles was considered as the true gradient. It was easy to notice that the LAGM gave a closer correlation with the true gradient than AGM and SPSA. For the calculated gradient correlation between FDM and the other three algorithms at the 10th iteration step, the SPSA marked by yellow line was 0.095, the AGM marked by gray line was 0.181, and the LAGM marked by blue line was 0.417. Similarly, at the 70th iteration step, the SPSA algorithm was 0.064, the AGM was 0.105, and the LAGM was 0.325. The results implied that the LAGM could effectively eliminate the pseudo correlation of the gradient. Fig. 9 was the convergence performance of the three algorithms within 70 iterations. In this example, LAGM and AGM yielded a final objective function value of about 7187. But LAGM showed a faster convergence speed, implying that LAGM had a higher computational efficiency for history matching process. The value of the final objective function obtained by SPSA was about 16601. Because SPSA required to perturb the variables synchronously for many times to obtain the approximate gradient, it often converged locally and required a higher computational cost. Also, it was tough to directly calculate the true gradient in history matching due to the huge numbers of reservoir parameters. Thus, parameterizing reservoir model and then modifying gradient might be a good choice for heavy reservoir history matching problem.
The final permeability estimates calculated by the three algorithms were presented in Fig. 10. The true model was obtained by FDM. By visual comparations, LAGM could reflect heterogeneous characteristics more accurate than AGM. In contrast, the permeability field obtained by SPSA gave more prominent heterogeneous characteristics. However, the obtained permeability field generated by SPSA was totally different as the true model in the upper left corner and lower right corner. The production predictions by the three algorithms were provided in Fig. 11. In this figure, the predictions of the initial model were marked by grey line. The observation data simulated by the true model was marked by the red circle. LAGM gave better history matching quality than the other two algorithms. However, for P3 well, no algorithms produced good results. Table 2 also showed that LAGM had higher calculation accuracy than the other two algorithms. This case mainly proved that the proposed method could effectively weaken the pseudo correlation and further improve the quality of history matching by using the single-well sensitive area and modified gradient.
3.2 Field Case
The Brugge oilfield was developed by TNO as a benchmark for closed-loop reservoir management test. Fig. 12 was the top structure of this reservoir. The grid represents the reservoir was divided into 139 × 48 × 9 grids, and the number of active grids is 44550. There are 20 production wells in the center and 10 injection wells at the edge to replenish energy within a lifetime of 7300 days. The injection rate was set to be 4000 m3/day and the production rate was 600 m3/day. In the history matching process, the variance of the measuring error was a fixed value of 0.05 over the observation data. The grid planar permeability was set to be the estimated parameter. Similar to the previously discussed example, 40 reservoir parameters were applied. Fig. 14a was the log-permeability field of the average model.
Then, FMM was introduced to calculate the flight time and to determine the single-well sensitive area. The sensitivity area of the first layer of the average model was shown in Fig. 13. The calculated flight time between grids with very long spatial distance was small due to the influence of strongly heterogeneous grids, resulting in noises. Therefore, we artificially eliminated some small noise in each area to improve the quality of history matching. Fig. 14 was the posterior estimate of permeability parameters generated by LAGM, AGM and SPSA algorithms at the 300th iteration step. LAGM and AGM gave very similar results in the MAP estimates, while SPSA gave more prominent heterogeneous characteristics. Fig. 15 was the matched production data of six wells. The production performances simulated by the average model diverted from the observations. LAGM algorithm showed better prediction results compared to other two algorithms. Table 3 presented the calculation accuracy of the three algorithms, indicating that LAGM had higher calculation accuracy in applications.
In this paper, we firstly reviewed the theorem of FMM and briefly introduced the sensitive area division method. Subsequently, we established a history matching mathematical model based on Bayesian framework. A parameterization method was used to reduce the dimensionality of original model parameters. By integrating the sensitive area information, the gradient of history matching could be better modified. The proposed approach was actually a two-step history matching method and no numerical simulation runs were required in the pretreatment process. To evaluate the feasibility of the presented method, a conceptual case of water drive in a heterogeneous two-dimensional reservoir and a real field case were tested. The results showed that the presented approach had higher accuracy and computational efficiency as compared to non-localized gradient method and gradient free algorithm. Moreover, the presented idea could also be used in other history matching algorithms, such as ENKF, ES-MDA, and so on. It was worth noting that parameterizing reservoir model and then modifying gradient may be a good choice for heavy reservoir history matching problem.
Funding Statement: This study was supported by National Natural Science Foundation of China (Nos. 52104017, 51874044, 51922007), Southern Marine Science and Engineering Guangdong Laboratory (Zhanjiang) (No. zjw-2019-04).
Conflicts of Interest: The authors declare that they have no conflicts of interest to report regarding the present study.
Appendix A: SPSA algorithm
As an approximate-gradient iterative optimization algorithm, SPSA can ensure that the optimization direction is always the uphill direction through the simultaneous disturbance of components of variable. For the lth iteration, the stochastic gradient of the target function in Eq. (11) can be calculated by
where is the random column vector that satisfies the Bernoulli distribution; is the positive coefficient controlling perturbation size with the fixed value of 0.101. The variable is updated through steepest descent
where αl+1 is the iteration step size with the fixed value of 0.5; the average gradient is defined by
where p is the number of perturbation at each iteration with a fixed value of 5.
|This work is licensed under a Creative Commons Attribution 4.0 International License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.|