Ensemble-based reservoir characterization method using multiple Kalman gains and dynamic data selection

Abstract

The present disclosure relates to an ensemble-based reservoir characterization method using multiple Kalman gains and dynamic data selection. The method includes preparing available data; generating initial ensembles by using the prepared static data; clustering and separating the generated initial models on the basis of a distance-based method; selecting the dynamic data; dynamically simulating the selected dynamic data by using the generated ensembles; calculating multiple Kalman gains by using initial models clustered in the same group as the selected dynamic data; updating ensemble members by means of the selected dynamic data and the multiple Kalman gains; and predicting a movement of a reservoir by using the updated models, and evaluating uncertainty thereof. Therefore, multiple Kalman gains are calculated and a final model is obtained using the selected dynamic data, and a reliable uncertainty evaluation and a future movement prediction can be performed within a short time by using the final model.

Claims

1. An ensemble-based reservoir characterization method using multiple Kalman gains and dynamic data selection, the method comprising: preparing available data including static and dynamic data; generating initial ensembles by using the prepared static data; clustering and separating the generated initial ensembles on the basis of a distance-based method; selecting data among all the dynamic data; dynamically simulating the selected dynamic data by using the generated initial ensembles; calculating multiple Kalman gains by using initial models clustered in the same group as the selected dynamic data; updating ensemble members by means of the selected dynamic data and the multiple Kalman gains; and predicting a movement of a reservoir by using the updated ensemble members, and evaluating uncertainty thereof and wherein at the updating the ensemble members by using an ensemble smoother, the state vector v, representing the ensembles, is expressed by the following equation: $y_{t_i,k,g}={\left[\begin{array}{c}{m}_{t_i}^s\\ m_{t_i}^d\\ d\end{array}\right]}_{k,g}$ where, g represents a group a model belongs to and is selected from 1 to the number of clusters (N.sub.c), k represents the model in the group g and is selected from 1 to the total number of models (N.sub.g) in the group g, thus yti,k,g represents a state vector of k-th model in the group g at an initial time ti, m.sub.t.sub.i.sup.a and m.sub.t.sub.i.sup.d represent static and dynamic parameters at the initial time, ti, respectively, and d represents dynamic data prediction; and the state vector is assimilated by using the multiple Kalman gains and the dynamic data selection by the following equation:
y.sub.k,g.sup.a=y.sub.k,g.sup.p+K.sub.g(d.sub.k,gHy.sub.k,g.sup.p) where, superscripts a and p indicate an assimilated and a previous states, respectively, y.sub.k,g.sup.a and y.sub.k,g.sup.p mean an assimilated and a previous state vectors of the k-th model in the group g, respectively, H represents a measurement operator matrix, and dk,g represents the dynamic data predication of the k-th model in the group g.

2. The method of claim 1, wherein at the calculating the multiple Kalman gains, the multiple Kalman gains are calculated by the following equation:
K.sub.g=C.sub.Y,g.sup.pH.sup.T(HC.sub.Y,g.sup.pH.sup.T+C.sub.D).sup.1 where, K.sub.g indicates a Kalman gain representing the group g among the multiple Kalman gains, C.sub.Y,g.sup.p indicates a previous estimate error covariance of the group g, H represents a measurement operator matrix, and C.sub.D represents an observation error covariance, and H.sup.T means a transpose of the measurement operator matrix.

3. The method of claim 2, wherein calculating the multiple Kalman gains includes calculating the multiple Kalman gains by group, after clustering the initial ensembles.

4. The method of claim 3, wherein the distance-based method is obtained by: obtaining a distance matrix; converting the distance matrix into coordinate points in the metric space; and classifying adjacent coordinate points into the same group.

5. The method of claim 4, wherein the obtaining the distance matrix is performed by using a Minkowski equation, wherein definition of a distance is performed for the static data or the selected dynamic data; the converting the distance matrix into the coordinate points in the metric space is performed by using multidimensional scaling; and the classifying the adjacent coordinate points into the same group is performed by using K-means algorithm.

6. The method of claim 1, wherein the selected dynamic data uses logically meaningful dynamic data only among all the available dynamic data, wherein the meaningful dynamic data have unique values different from the dynamic data.

Description

DESCRIPTION OF DRAWINGS

(1) FIG. 1 is a flow chart illustrating a process for prediction of reservoir characterization by using the standard ES.

(2) FIG. 2 is a flow chart illustrating a process for the ES with multiple Kalman gains and dynamic data selection according to an embodiment of the present disclosure.

(3) FIGS. 3a and 3b are diagrams illustrating the workflow of the ES and EnKF, respectively, with observed dynamic data (FIG. 3c) according to an embodiment of the present disclosure, wherein FIG. 3a is the ES; FIG. 3b is the EnKF; and FIG. 3c is the observed dynamic data.

(4) FIG. 4 is a diagram illustrating a procedure of clustering static models and obtaining the Kalman gain for each group (the concept of multiple Kalman gains) according to an embodiment of the present disclosure.

(5) FIG. 5 is a diagram of top structure of production wells illustrating PUNQ-S3 model expressed on the X-Y plane according to an embodiment of the present disclosure.

(6) FIGS. 6a, 6b, 6c, and 6d are diagrams illustrating horizontal permeability (static model) according to an embodiment of the present disclosure, wherein FIG. 6a illustrates the reference field regarded as a right answer; FIG. 6b illustrates an average of the initial models; FIG. 6c illustrates an average of results of the standard ES; and FIG. 6d illustrates an average of results of the ES with multiple Kalman gains.

(7) FIG. 7 is a diagram illustrating the similarity of the three static models in the same cluster for three clusters as a result of the distance-based clustering according to an embodiment of the present disclosure.

(8) FIGS. 8a, 8b, 8c, and 8d illustrate the histograms of horizontal permeability from the first layer of the PUNQ-S3 model according to an embodiment of the present disclosure, wherein FIG. 8a is a graph of the reference field; FIG. 8b is a graph of the initial model; FIG. 8c is a graph of the standard ES; and FIG. 8d is a graph of the ES with multiple Kalman gains.

(9) FIGS. 9a, 9b, and 9c illustrate prediction of each reservoir performance from the updated models of the standard ES (the first row) and the ES with multiple Kalman gains (the second row) according to an embodiment of the present disclosure, wherein FIG. 9a is a graph illustrating cumulative oil production (COP); FIG. 9b is a graph illustrating cumulative gas production (CGP); and FIG. 9C is a graph illustrating cumulative water production (CWP).

(10) FIGS. 10a, 10b, 10c, and 10d illustrate horizontal permeability of all five layers of the PUNQ-S3 model according to an embodiment of the present disclosure, wherein FIG. 10a is a diagram illustrating the reference field; FIG. 10b is a diagram illustrating averages of initial models; FIG. 10c is a diagram illustrating averages of the ES with selected dynamic data; and FIG. 10d is a diagram illustrating averages of the ES with multiple Kalman gains and selected dynamic data.

(11) FIGS. 11a, 11b, 11c, and 11d illustrate the histograms of horizontal permeability of a PUNQ-S3 model according to an embodiment of the present disclosure, wherein FIG. 11a is a diagram of the reference field; FIG. 11b is a diagram of the initial model; FIG. 11c is a diagram of the ES with selected dynamic data; and FIG. 11d is a diagram of the ES with multiple Kalman gains and selected dynamic data.

(12) FIGS. 12a and 12b illustrate predictions of CWP from the updated models according to an embodiment of the present disclosure, wherein FIG. 12a is a graph for the ES with selected dynamic data; and FIG. 12b is a graph for the method of the present disclosure, ES with multiple Kalman gains and selected dynamic data.

MODE FOR INVENTION

(13) Hereinafter, exemplary embodiments of the present disclosure will be described in detail with reference to the accompanying drawings.

(14) FIG. 2 is a flow chart illustrating a process for the ES with multiple Kalman gains and dynamic data selection according to an embodiment of the present disclosure.

(15) In step S202, available data are prepared. The available data include static and dynamic data.

(16) In step S204, initial ensembles are generated by using the static data.

(17) In step S206, the initial ensembles are classified by clustering based on a distance-based method, wherein the initial ensembles are generated in step S204. The initial ensembles similar to each other are bound into a same group through this.

(18) In step S208, meaningful data are selected among the dynamic data prepared in step S202.

(19) In step S210, dynamic simulation for the dynamic data selected in step S208 is performed by using the initial ensembles generated in step S204.

(20) In step S212, multiple Kalman gains along with selected dynamic data in step S208 are calculated by using the initial ensembles bound into the same group in step S206.

(21) In step S214, ensemble members are updated with selected dynamic data in step S208 and multiple Kalman gains in step S212.

(22) In step S216, reservoir performances are predicted or uncertainty is quantified using the updated ensemble members.

(23) The steps for calculating multiple Kalman gains and updating a static model by using dynamic data will be described more in detail with reference to the following drawings.

(24) FIGS. 3a and 3b are diagrams illustrating the workflow of the ES and EnKF, respectively, with observed dynamic data (FIG. 3c) according to an embodiment of the present disclosure, wherein FIG. 3a is the ES; FIG. 3b is the EnKF; and FIG. 3c is the observed dynamic data.

(25) Referring to FIGS. 3a, 3b, and 3c, the ES assimilates all available dynamic data at a time regardless of observed time, which is called a global update (FIG. 3a). Accordingly, only one forecast step is required from an initial time, t.sub.i, to a last observed time, t.sub.i, as shown in equation 1.

(26) $\begin{array}{cc}{y}_{t_i,k,g}={\left[\begin{array}{c}{m}_{t_i}^s\\ m_{t_i}^d\\ d\end{array}\right]}_{k,g}& \left[\mathrm{Equation}2\right]\end{array}$

(27) where, m.sub.t.sub.i.sup.d and m.sub.t.sub.i.sup.d represent dynamic variables at the initial time, t.sub.i, and last time, t.sub.i, respectively, m.sub.t.sub.i.sup.x represents a static variable at the initial time, t.sub.i, and d represents dynamic data prediction from initial model to be compared with observed dynamic data. In the case of the EnKF, a recursive update is used, and each model is simulated until the next observed time instead of the last observed time (FIG. 3b). Other equations and assumptions except the update concept for assimilation are the same with those of the ES.

(28) After the forecast step, each ensemble is expressed in the state vector as equation 2 to perform assimilation.

(29) $\begin{array}{cc}\left\{\begin{array}{c}{m}_{t_i}^d\\ d\end{array}\right\}=f\left(m_{t_i}^s,m_{t_i}^d\right).& \left[\mathrm{Equation}1\right]\end{array}$

(30) where, g, ranging from 1 to the number of total clusters (N.sub.c), represents the group the model belongs to, k represents the models, which belong to the group g ranging from 1 to the number of total models in the group g (N.sub.g), and y.sub.ti,k,g represents the state vector of k-th model in the group g at initial time, t.sub.i. As shown in equation 2, the state vector of the k-th ensemble member that belongs to the group g is composed of three elements (m.sub.ti.sup.x, m.sub.ti.sup.d, d).

(31) In addition, the Kalman gain that minimizes the estimate error covariance is calculated as equation 3.
K.sub.g=C.sub.Y,g.sup.pH.sup.T(HC.sub.Y,g.sup.pH.sup.T+C.sub.D).sup.1[Equation 3]

(32) where, K.sub.g indicates the Kalman gain representing the group g among multiple Kalman gains, C.sub.Y,g.sup.p indicates the previous estimate error covariance of the group g, H and C.sub.D represent the measurement operator matrix and observation error covariance, respectively, and H.sup.T means a transpose of the measurement operator matrix.

(33) The Kalman gain has a similar role with the sensitivity matrix in gradient-based optimization methods as shown in equation 4.
m.sub.i.sup.u=m.sub.i.sup.uc+C.sub.MG.sup.T(GC.sub.MG.sup.T+C.sub.D).sup.1(d.sub.ig(m.sub.i.sup.uc))[Equation 4]

(34) where m.sub.i.sup.u and m.sub.i.sup.uc represent the assimilated and previous i-th model vectors, respectively, C.sub.M means the covariance matrix of the static parameters, and G means the sensitivity matrix.

(35) As the sensitivity matrix has different values depending on model parameters, the Kalman gain should be also calculated depending on the models as well. However, previous researches on the ensemble-based method used a unified single Kalman gain for all ensemble members.

(36) In the present disclosure, the Kalman gain for a corresponding group, K.sub.g, is calculated by clustering similar ensembles as shown in equation 3 above. Equation 5 below is exactly same with one of the standard ES except multiple Kalman gains.

(37) It is used for equation 5 to assimilate the state vector y.sub.ti,k,g in equation 2 above.
y.sub.k,g.sup.a=y.sub.k,g.sup.p+K.sub.g(d.sub.k,gHy.sub.k,g.sup.p)[Equation 5]

(38) where, superscripts a and p indicate an assimilated and a previous states, respectively, y.sub.k,g.sup.a and y.sub.k,g.sup.p mean an assimilated and a previous state vectors of the k-th model in the group g, respectively, and d.sub.k,g represents the dynamic data prediction of the k-th model in the group g.

(39) In the ensemble-based method, it is assumed that the mean of ensembles is true value and the static parameters follow Gaussian distribution. Therefore, the ensemble-based method is difficult to apply to improper initial ensemble models due to lack of geological information and to non-Gaussian fields such as channelized reservoirs.

(40) The proposed method in the present disclosure improves the assimilated results by using multiple Kalman gains and selected dynamic data in equation 5. The concept of selected measurement data is also used in equation 3 above for the calculation of multiple Kalman gains. Since selective use of measurement data is a simple and intuitive concept, it can be one of solutions for reliable history matching. The ensemble-based method can assimilate various types of dynamic data such as pressure, production, and injection data. However, some data are not sensitive for model parameters and give confusing information for assimilation. Therefore, only meaningful production data are selected logically instead of all the available dynamic data.

(41) These two concepts (multiple Kalman gains and selected dynamic data) are coupled with the standard ES and it is verified for PUNQ-S3 model.

(42) FIG. 4 is a diagram illustrating a procedure of clustering static models and obtaining the Kalman gain for each group (the concept of multiple Kalman gains) according to an embodiment of the present disclosure.

(43) Referring to FIG. 4, the Kalman gain for each cluster can be calculated by clustering initial models with the distance-based method. After initial models are generated by integration of static data, a distance-based method is utilized to classify the initial models based on dissimilarity. The process for obtaining clustered covariance typically follows three steps: defining a distance, converting the models into the metric space using the defined distance, and clustering algorithm for grouping the models. At first, a Minkowski equation is utilized to obtain the distance matrix, which is defined as dissimilarity between two realizations, as shown in equation 6.

(44) $\begin{array}{cc}{d}_{\mathrm{ab}}=\sqrt[p]{\underset{i=1}{\overset{N}{.Math.}}{.Math.x_{\mathrm{ai}}-x_{\mathrm{bi}}.Math.}^p}& \left[\mathrm{Equation}6\right]\end{array}$

(45) where, d.sub.ab represents the distance between models a and b, x represents a parameter compared for the distance calculation for each model and can be anything as long as it is related to a characterization of each model, x.sub.ai and x.sub.bi represent the values of i-th x parameter of the initial models a and b, respectively, and |x.sub.aix.sub.bi|.sup.p represents the absolute distance (positive value) between the models a and b with respect to the x parameter.

(46) As described above, some researchers defined the x, which is a parameter for comparison of models, by using static data such as rock facies or permeability. Others employed dynamic data such as OPR or water breakthrough time for x. Although all of them successfully distinguished reservoir models, distance calculation from dynamic data requires additional simulation time. Clustering time can be reduced with the dissimilarity calculation defined by using static data according to a present embodiment.

(47) Then, the distance matrix is converted into coordinates in n-dimensional (n-D) metric space by MDS, which is one of popular dimensional reduction methods. It depicts the models in the metric space while maximally conserving the original distance among models. The points in the metric space in FIG. 4 stand for each reservoir model and the more different the models are, the more increased the distances in the metric system are.

(48) Finally, K-mean algorithm is applied to classify neighbor points into the same group. This classification consists of two steps of: selecting centroids and allocating the points to the closest centroid. The result of clustering is sensitive for location of the first set of centroids and the number of cluster. Silhouette index and sensitivity analysis are usable algorithms for the determination of the number of clusters.

(49) FIG. 5 is a diagram of top structure of production wells illustrating PUNQ-S3 model expressed on the X-Y plane according to an embodiment of the present disclosure.

(50) Referring to FIG. 5, the proposed method of the present disclosure, the ES with multiple Kalman gains and selective dynamic data, is applied to a widely known benchmark problem, PUNQ-S3 model, to confirm the effect thereof.

(51) The grid of the PUNQ-S3 model consists of 19285 with x=y=180 m and is based on the corner-point system as illustrated in FIG. 5. In addition, it is also surrounded by a fault expressed in thick black line on the southeast part and bounded by an aquifer expressed in thick line on the northwest, which affects the number of active grids in each layer. The numbers of active grids are 396 on the first and second layers and 323 on the other layers, respectively. In FIG. 5, there are six production wells which have static data and four kinds of dynamic data: well OPR (WOPR), well bottom-hole pressure (WBHP), well gas-oil ratio (WGOR), and well WCT (WWCT).

(52) In the experiment of the present disclosure, the total simulation period for the PUNQ-S3 model is 16.5 years, which is separated into the first 8 years of assimilation and the following 8.5 years of prediction, wherein current point is eighth year. During the first 8 years, well tests were carried out during the first year and then wells were shut-in for 3 years, and for the next 4 years, production was achieved. ECLIPSE 100 of Schlumberger is used for dynamic reservoir simulation.

(53) FIGS. 6a, 6b, 6c, and 6d are diagrams illustrating horizontal permeability (static model) according to an embodiment of the present disclosure, wherein FIG. 6a illustrates the reference field regarded as a right answer; FIG. 6b illustrates an average of the initial models; FIG. 6c illustrates an average of results of the standard ES; and FIG. 6d illustrates an average of results of the ES with multiple Kalman gains.

(54) Referring to FIGS. 6a, 6b, 6c, and 6d, FIG. 6a is the reference field according to an embodiment of the present disclosure as the horizontal permeability in the first layer of the PUNQ-S3 model; FIG. 6b is the average of 100 initial ensemble models generated with geostatistical method from the static data of the six production wells of FIG. 5; FIG. 6c shows the result of the standard ES; and FIG. 6d shows the result of the ES with multiple Kalman gains. Meanwhile, the porosities for all layers are generated by sequential Gaussian simulation (SGS) from static geological information available, then initial horizontal permeability models are estimated with the porosity models by co-SGS. The correlation coefficient for porosity and horizontal permeability used in the present disclosure is 0.83, while geostatistics are implemented by SGeMS (Stanford geostatistical modeling software).

(55) Prior to applying the proposed method of the present disclosure, the effect of multiple Kalman gains was investigated first. Here, only horizontal permeability from the first layer is used for model parameter as shown in FIG. 6a. The distance for clustering is defined by WGOR at PRO-12 well which has a wide response range for initial ensembles.

(56) FIG. 7 is a diagram illustrating the similarity of the three static models in the same cluster for three clusters as a result of the distance-based clustering according to an embodiment of the present disclosure. Referring to FIG. 7, each point in FIG. 7 represents one of initial models and initial models in the same cluster have similar distribution of static parameters. Multiple Kalman gains are calculated by using the result of clustering and are used in assimilation. The updated models with the standard ES from 100 initial ensembles show an obvious overshooting problem as illustrated in FIG. 6c. However, as illustrated in FIG. 6d, the ES with multiple Kalman gains mitigates such the problem and conserves the major pattern of the reference field without the concept of selective use of dynamic data.

(57) FIGS. 8a, 8b, 8c, and 8d illustrate the histograms of horizontal permeability from the first layer of the PUNQ-S3 model according to an embodiment of the present disclosure, wherein FIG. 8a is a graph of the reference field; FIG. 8b is a graph of the initial model; FIG. 8c is a graph of the standard ES; and FIG. 8d is a graph of the ES with multiple Kalman gains.

(58) As illustrated in FIG. 8a, the histogram of the reference field has a highly negatively skewed distribution. In FIG. 8b, the histogram of the initial ensemble model has normal distribution due to the SGS algorithm. In FIG. 8c, the histogram of the standard ES shows the obvious overshooting problem and maximum and minimum values about 36.5 and 17.3 for the logarithm of permeability, respectively, which are unrealistic values for permeability. As a result, the standard ES is unable to give a reliable reservoir characterization for the PUNQ-S3. However, the histogram of the ES with multiple Kalman gains in FIG. 8d solves the overshooting problem and has a relative frequency nearby 7 log-millidarcy (log-mD). Consequently, this shows that the ES with multiple Kalman gains has an effect of improving the ES for assimilation.

(59) FIGS. 9a, 9b, and 9c illustrate prediction of each reservoir performance from the updated models of the standard ES (the first row) and the ES with multiple Kalman gains (the second row) according to an embodiment of the present disclosure, wherein FIG. 9a is a graph illustrating COP; FIG. 9b is a graph illustrating CGP; and FIG. 9C is a graph illustrating CWP.

(60) Referring to FIGS. 9a, 9b, and 9c, the thin solid lines represent the result of each updated model and the thick solid line stands for the reference field's production. The dot lines from the top are P10, P50, and P90 of the thin solid lines, respectively.

(61) Since the ES is assimilated with wrong models, as is confirmed in FIGS. 6c and 8c, the results from the standard ES show unreliable future reservoir performances accordingly, as illustrated in the first row of FIGS. 9a, 9b, and 9c. The thin solid lines for COP, CGP, and CWP are far off from the true values of thick solid line. However, the productions from the ES with multiple Kalman gains give much more improved and reliable ranges including true values (the second row of FIGS. 9a, 9b, and 9c). This is due to use of multiple Kalman gains based on the clustering instead of the inaccurate one Kalman gain.

(62) FIGS. 10a, 10b, 10c, and 10d illustrate horizontal permeability of all five layers of the PUNQ-S3 model according to an embodiment of the present disclosure, wherein FIG. 10a is a diagram illustrating the reference field; FIG. 10b is a diagram illustrating averages of initial models; FIG. 10c is a diagram illustrating averages of the ES with selected dynamic data; and FIG. 10d is a diagram illustrating averages of the ES with multiple Kalman gains and selected dynamic data.

(63) The proposed method of the present disclosure, the ES with multiple Kalman gains and selected dynamic data, is applied to the 3D reference fields like FIG. 10a. For multiple Kalman gains, horizontal permeability, which is the static model parameter itself, is employed to the distance calculation. For selected dynamic data concept, Table 1 shows the selected dynamic data logically among all observed data and observation noises. For example, shut-in BHP data among WBHP are eliminated from the data for assimilation because they are very similar for each initial model. Otherwise, they may increase calculation time in assimilation by enlarging matrix size and decrease assimilation reliability by confused information.

(64) TABLE-US-00001 TABLE 1 Time WBHP WGOR WWCT (Days) (number) (number) (number) 274 6 1,642 6 1,826 6 6 2,008 6 2,192 6 6 6 2,373 6 2,557 6 6 6 2,572 6 2,738 6 2,992 6 6 6 Total 24 36 42 Observation 3 bar 10% for less 0.01 noise than 90 m.sup.3/m.sup.3, 25% for more than 90 m.sup.3/m.sup.3

(65) When the average of updated horizontal permeability for initial models is checked, no overshooting problem is found in cases of the ES with selective data and the proposed method (the ES with multiple Kalman gains and selective data) as shown in FIGS. 10c and 10d. In addition, the averages of updated models show proper specific patterns of the reference fields. The previous studies of a standard EnKF for the PUNQ-S3 or the studies of the standard ES above (see FIGS. 6c and 8c) showed overshooting problem, but the problem is resolved with selective use of the dynamic data for the ES. Thus, it is confirmed that selective use of the dynamic data for the ensemble-based method is important and effective.

(66) FIGS. 11a, 11b, 11c, and 11d illustrate the histograms of horizontal permeability of the PUNQ-S3 model according to an embodiment of the present disclosure, wherein FIG. 11a is a diagram of the reference field; FIG. 11b is a diagram of the initial model; FIG. 11c is a diagram of the ES with selected dynamic data; and FIG. 11d is a diagram of the ES with multiple Kalman gains and selected dynamic data.

(67) The histogram of the reference field has a high relative frequency around 7 of log-mD as illustrated in FIG. 11a. In the case of the ES with selective dynamic data, overshooting values in the histogram of the standard ES as illustrated in FIG. 8c disappear as illustrated in FIG. 11c due to the effect of selective use of observed data. In addition, the histogram of the proposed method of the present disclosure, the ES with multiple Kalman gains and selected dynamic data, reveals no overshooting problem as shown in FIG. 11d. Furthermore, it describes properly the distribution of the reference field having a high relative frequency around 7 of log-mD.

(68) FIGS. 12a and 12b illustrate predictions of CWP from the updated models according to an embodiment of the present disclosure, wherein FIG. 12a is a graph for the ES with selected dynamic data; and FIG. 12b is a graph for the method of the present disclosure.

(69) The result from the ES with selected data as illustrated in FIG. 12a gives better prediction than those from the standard ES as illustrated in first row of FIGS. 9a, 9b, and 9c. However, uncertainty thereof is so highly biased that the true value (thick solid line) is positioned out of the uncertainty range between P10 and P90 as illustrated in FIG. 12a. The results from the method of the present disclosure, the ES with multiple Kalman gains and selected dynamic data, show reliable uncertainty quantification that the actual production is positioned in the middle of uncertainty range. Consequently, reasonable decision making and development plan can be established through prediction of future production from the proposed method only.

(70) Meanings according to an embodiment of the present disclosure are as follows. Firstly, the method of the present disclosure is tested to one of the most famous benchmark models, PUNQ-S3. The benchmark model is known to be difficult to characterize even for EnKF due to the overshooting problem. It is very hard to apply the ES for PUNQ-S3, since the ES, despite of the faster calculation speed thereof, is more vulnerable to the overshooting problem than the EnKF. Actually, the standard ES shows serious overshooting problem even in 2D case as can be seen in FIGS. 6c and 8c. However, the method of the present disclosure carried out successfully the characterization of PUNQ-S3 by using the two concepts of multiple Kalman gains and selected observed dynamic data.

(71) Secondly, initial models are successfully grouped for 3D model. Initial models are classified by a Minkowski distance, MDS, and K-means clustering using a permeability value itself, or selected dynamic data. Reliable clustering is prerequisite for a proper calculation of multiple Kalman gains.

(72) Thirdly, the ES with multiple Kalman gains without introducing selective dynamic data concept resolves the overshooting problem of the standard ES and conserves direction of high permeability in the reference field. The concept of multiple Kalman gains can increase reliability of assimilation because a proper Kalman gain can be coupled with each ensemble instead of only one Kalman gain.

(73) Fourthly, selective use of the observed dynamic data mitigates the overshooting problem that the standard ES has. However, the ES with selective dynamic data has some limitation that the reliability of future production from the updated models is low.

(74) Fifthly, the ES with multiple Kalman gains and selective dynamic data (the method of the present disclosure) resolved the problems arising from the standard ES or the ES with selective dynamic data concept only. Future production from the updated fields gives an appropriate uncertainty range including the true production. Thus, the present disclosure can expand the applicability of the ES for complex 3D benchmark model and can be utilized for reasonable decision making.

(75) Sixthly, the method of the present disclosure as described above can be processed by a general computing apparatus, wherein the computing apparatus can be materialized by using any appropriate combination of at least one of processing circuit executing a software module stored in a memory such as a microprocessor, a digital signal processor (DSP), a microcontroller, a portable wireless phone, a personal digital assistant (PDA), and a paging device; or a firmware in which a software and/or an application stored in hardware and the memory is maintained. For example, the specifications of the computer used in the present research are 3.4 GHz CPU consisting of eight cores and 32 GB RAM consisting of four 8 GB memory modules. In addition, all calculations were done in MATLAB and ECLIPSE 100 is also implemented through MATLAB.

(76) Seventhly, the method of the present disclosure not only improves the assimilation reliability but also reduces sharply time required for history matching. In comparison of time required for history matching, a total of 5,687 minutes are necessary for the case of the EnKF due to 45 iterative assimilations while a total of 128 minutes is necessary for the proposed ES according to the present disclosure since one assimilation is only required. Consequently, the proposed ES according to the present disclosure saves time by 97.7% compared with the EnKF. The preprocessing stage of the present disclosure needs only additional four minutes because distance-based clustering is not an inverse modeling algorithm. Therefore, the proposed method is almost similar with the standard ES in terms of simulation cost.

(77) Although the content of the present disclosure above has been described with reference to embodiments shown in drawings for illustrative purposes, those skilled in the art will appreciate that various modifications and equivalent other embodiments are possible. Therefore, the real technical protection scope of the present disclosure should be defined by the technical idea as disclosed in the accompanying claims.