Integration of an Iterative Update of Sparse Geologic Dictionaries with ES-MDA for History Matching of Channelized Reservoirs

Severe Storm Research Center, Ewha Womans University, 52 Ewhayeodae-gil, Seodaemun-gu, Seoul 03760, Republic of Korea Department of Climate and Energy Systems Engineering, Ewha Womans University, 52 Ewhayeodae-gil, Seodaemun-gu, Seoul 03760, Republic of Korea Petroleum and Marine Research Division, Korea Institute of Geoscience and Mineral Resources, 124 Gwahak-ro, Gajeong-dong, Yuseong-gu, Daejeon 34132, Republic of Korea Department of Energy Resources Engineering, Seoul National University, 599 Gwanak-gu, Seoul 151-744, Republic of Korea


Introduction
Calibration of a subsurface system is an essential process to forecast fluid behaviors in a variety of geoenvironments such as aquifers, geothermal reservoirs, and petroleum reservoirs.History matching is an inverse process to find reservoir model parameters honoring observations by integration of static (e.g., core, logging, and seismic) and dynamic data (e.g., oil and gas rate, water cut, bottomhole pressure, and subsidence/uplift) [1].Ensemble-based data assimilation approaches have been successfully utilized for history matching to provide subsurface models that are well conditioned to observations.For example, the ensemble Kalman filter (EnKF) [2][3][4][5], ensemble smoother (ES) [6,7], and ES with multiple data assimilation (ES-MDA) [8][9][10].However, the ensemble-based data assimilation approaches have difficulty in preserving non-Gaussian distributions of model parameters such as lithofacies [11][12][13][14].In the ensemble-based data assimilation approaches, model parameters lose the non-Gaussianity of their original distributions that are initially constrained and the distributions of the model parameters get close to Gaussian ones.Shin et al. [14] and Zhou et al. [15] suggested using normal score transform in the ensemble-based data assimilation approaches to preserve non-Gaussian distributions of model parameters.Non-Gaussian model parameters are transformed into Gaussian model parameters using normal score transform, and then finally, updated model parameters are backtransformed.Moreover, transformation can take advantage of parameterization if the number of essential transformed parameters is smaller than the number of original parameters in terms of saving computational cost and figuring out main features of parameters.For example, the discrete cosine transform (DCT) [16][17][18][19][20][21], fast Fourier transform [22,23], grid connectivity transform [24], level set [13,25], and sparse geologic dictionaries [26,27] have been applied to reservoir characterization.In particular, Fourier transform-based methods such as DCT are capable of capturing essential traits such as main shapes and patterns of a facies channel reservoir [16,17] but reveal a deficiency in describing a crisp contrast among different facies because of data loss from inverse transformation [28].
Sparse coding refers to the process of computing representation coefficients based on the given signal and dictionaries [29].In sparse coding, the dictionaries indicate groups of features capable of brief expressions to represent various phenomena in the environment [30].In geological modeling, sparse geologic dictionaries are used to represent models with sparse linear combinations of basis vectors that are essential geologic features of a reservoir [29].Extracting essential geologic features and reducing the number of reservoir variables can be accomplished using sparse coding, thereby facilitating ensemble-based history matching [27].Aharon et al. [29] showed the efficacy of K-singular value decomposition (K-SVD) resulting from the accelerated convergence for image reconstruction, which led Sana et al. [31] to build an archive of essential geologic features called sparse geologic dictionaries from thousands of static reservoir models using K-SVD and calibrate reservoir models with the dictionaries using EnKF.One drawback of K-SVD is its large size of sparse geologic dictionaries.References showed that sparse coding with a transformation of parameter space could reduce both computational complexity and costs that are simultaneously required for model calibration [26,27,32].In this study, we note that the previous works have not considered the quality of sparse geologic dictionaries.The quality indicates how well reservoir models can be properly reconstructed by prototypes of dictionaries.Also, we expect an improvement in the history-matching performance by enhancing the quality of sparse geologic dictionaries.
This study proposes a hybridized ES-MDA algorithm that implements sparse coding in a transformed space to outperform previous history-matching methods by providing more accurate reconstructions of highly non-Gaussian model parameters.The proposed algorithm transforms multimodal facies into coefficients of discrete cosine functions using DCT.Invoking DCT is followed by iterating K-SVD for updating sparse geologic dictionaries.In each assimilation of ES-MDA, the combination of DCT and iterative K-SVD is performed to update the dictionaries and improve the quality of reservoir models.For brevity, the proposed algorithm with updated sparse geologic dictionaries is called ES-MDA-DCT-i-K-SVD in this paper.The performance of ES-MDA-DCT-i-K-SVD is tested with applications for channelized gas reservoirs and is compared with those of four ES-MDA algorithms: conventional ES-MDA, ES-MDA coupled with DCT (called ES-MDA-DCT in this paper), ES-MDA coupled with K-SVD (called ES-MDA-K-SVD in this paper), and ES-MDA coupled with DCT and K-SVD (called ES-MDA-DCT-K-SVD in this paper).

Methodology
The novelty of the proposed algorithm ES-MDA-DCT-i-K-SVD is the integration of ES-MDA (Section 2.1), the dimensionality reduction of the parameter space using DCT (Section 2.2), and construction (Section 2.3) and updating (Section 2.4) of geologic dictionaries using sparse coding in the reduced space.Section 2.5 describes how the methods operate in the framework of ES-MDA-DCT-i-K-SVD in a complementary manner.

ES-MDA.
The goal of history matching can be formulated as where J is the objective function of history matching and m is the state vector composed of reservoir variables (e.g., permeability and facies).
The typical form of J m for ensemble-based history matching is presented as [33] where m b is the state vector before update and the superscript b refers to background; B is the covariance matrix of m b ; d obs is the observed responses; d = f m is the dynamic vector composed of simulated responses obtained by running a reservoir simulator f for the state vector m; and R is the covariance matrix of observation error.Note that the right-hand side of ( 2) is the addition of background and observation error terms [33].Because m can contain any unknown variables such as facies indexes, coefficients of discrete cosine functions or dictionary coefficients depending on the type of algorithms were used in this study.∂J m /∂m = 0 can be used to derive the update equation for m as [8,33] where the subscript i refers to the ith ensemble member; C md is the cross-covariance matrix of m and d; C dd is the autocovariance matrix of d; α p is the coefficient to inflate C D , which is the covariance matrix of the observed data measurement error [8]; d unc is the observation data perturbed by the inflated observed data measurement error; and N ens is the ensemble size (i.e., number of reservoir models in the ensemble).Conventionally, ensemble-based history matching updates N ens reservoir models simultaneously.In (3), C md C dd + α p C D −1 refers to Kalman gain K, which is computed with regularization by SVD using 99.9% of the total energy in singular values [8].

Geofluids
The main difference between ES and ES-MDA is the update process of the state vector m.ES updates the state vector of each ensemble member using observation data measured at all time steps (Emerick and Reynolds, 2011).Compared to the single assimilation of ES, ES-MDA assimilates every state vector N a times using an inflated covariance matrix of measurement error [8,9].Here, N a is the number of assimilations in ES-MDA.
Definitions of C md and C dd are as follows: where m is the mean of state vectors and d is the mean of dynamic vectors.
In ES-MDA, α p is constrained to In ES, N a = 1 and α p=1 = 1 due to its single assimilation.The perturbed observation data d unc shown in ( 3) is computed as The second term on the right-hand side of ( 6) is the perturbation term, which reflects the uncertainty associated with

Extraction of Geologic Features Using Discrete Cosine
Transform.Discrete cosine transform (DCT) has been utilized as an image-processing tool for characterization of channelized reservoirs due to the periodicity of cosine functions [34].DCT converts parameters into coefficients of discrete cosine functions.The coefficients are sorted in descending order from the top left, capturing the overall trend of channel patterns, to bottom right, delineating details in channel patterns.Previous studies have shown that non-Gaussian channel patterns can be reproduced sufficiently via inverse transform of essential DCT coefficients [18,28,35].Updating the truncated DCT coefficients can yield a calibrated model set.Another advantage of DCT is the improvement in computational efficiency resulting from data compression, which is effective in constructing sparse geologic dictionaries described in Section 2.3.
Figure 1 illustrates how to extract geologic features from an image of a target channelized reservoir using a truncated DCT and reproduce the target reservoir through an inverse DCT (IDCT).Two images in the first row represent the physical state of sand and shale facies in the target reservoir, that is, the original image on the left and the reproduced image on the right.Let N grid and N DCT denote the number of gridblocks of the reservoir model and the number of essential DCT coefficients, respectively.Applying DCT to the original  1, filtering the coefficient state selects 465 components in the dotted triangle as essential ones.It seems that this small number of components is sufficient to restore the original image of the physical state (i.e., channel patterns) when comparing the two subfigures in the first row.

Construction of Geologic Dictionaries Using Sparse
Coding. Figure 2 presents a procedure of sparse coding to construct geologic dictionaries using static data in the DCT domain.Sparse coding starts from building a library matrix Y, which is a N para by N lib matrix.N para is the number of parameters in each reservoir model, and N lib is the number of reservoir models in Y. Equiprobable reservoir models are generated using a geostatistical technique.Figure 2(a) shows a training image used for creating initial channelized reservoir models (see Figure 2(b)) by invoking single normal equation simulation (SNESim) [36].In Figures 2(a That is, the data compression ratio is N grid /N DCT .In Figure 2, N grid /N DCT = 75 × 75/465 ≈ 12 1.Data compression using DCT reduces the dimension of the parameter space, thereby saving computing time required for sparse coding [37].Meanwhile, a sufficiently large N lib needs to be chosen to cover a variety of geologically plausible scenarios in Y. Previous investigators adopted N lib in the range of 1000 to 2000 [31,37].In this paper, N lib is a constant of 3000 for maintaining the diversity of library models.Note that generating N lib models is computationally inexpensive because all the models are static.No dynamic reservoir simulation is required for any static model. Library matrix Y is numerically decomposed into dictionary matrix D and weight matrix X using K-SVD: Y ≅ DX.D and X are an N para by N dict matrix and an N dict by N lib matrix, respectively, where N dict is the number of reservoir models in D. Thus, both Y and D are sets of reservoir models.Each column vector of D, called a realization in this paper, represents either an original or an encoded reservoir model.We expect that every realization exhibits distinguishable geological features in a well-organized dictionary.N dict is to be predetermined considering computational costs associated with the sparse coding process.In this paper, N dict is fixed as one third of N lib .As a rule of thumb,  where • F is the Frobenius norm and Y′ = DX.More specifically, matrix decomposition is performed by iterating sparse coding that is an orthogonal matching pursuit (OMP) [38,39] followed by K-SVD [29], as described by Sana et al. [31].The first step of matrix decomposition is to initialize D and X.The second step is to compute X with a fixed D using OMP.The third step is update D and X simultaneously using K-SVD through the following three equations.
where d r i is the ith column vector (i.e., realization) in D and x i is the ith row vector in X.
The right-hand side of ( 8) is rearranged as follows: where the subscript j ∈ 1, … , N dict indicates the pivot and E j is the discrepancy term.
To achieve ( 9), optimal d r j and x j are explored from j = 1 to j = N dict .The right-hand side of ( 9) is rearranged as where Ω j is the matrix of which every element is either zero or one.At each j, x j multiplied by Ω j returns a row vector of nonzero elements in x j .The size of Ω j is N lib by the number of nonzero elements in x j .E j Ω j is decomposed using SVD for determining d r j and x j at each j, which results in U ΔV T , where U is the first column vector from the leftsingular vectors' matrix, Δ is the first element of the singular values' matrix, and V is the first column vector from the right-singular vectors' matrix of E j Ω j , respectively.As a consequence, d r j is the first column vector of U and x j Ω j is the first column vector of V multiplied by the first diagonal element of Δ. x j is obtained by multiplying the inverse Ω j .Note that the updated d r j and x j are used for the calculation of (9) at the subsequent js.Performing (10) at all js completes the update of D and X.
After the update of the dictionary matrices, the second step of matrix decomposition is reinvoked to tune X for further achieving (7).Then, the tuned X and D are reimported to (8) for conducting the third step again.This sequence of the second and third steps is iterated until a convergence criterion is satisfied.The criterion could be sparsity of X, a threshold to accept the discrepancy shown in (7) or the maximum number of iterations.In this study, the criterion is the maximum number of iterations (set to ten) for all experiments.
In summary, Figures 2(a)-2(d) describe how to construct the original sparse geologic dictionary matrices Y, D, and X using static data.designation.Note that the diffusion is filtered out using a cutoff method in this study.We employ the arithmetic mean of the two facies indexes (i.e., 0.5) as the threshold to determine either shale (if facies index < 0.5) or sand (if facies index ≥ 0.5) in each gridblock.

Update of Sparse Geologic Dictionaries during Multiple
Data Assimilation.Once original sparse geologic dictionary matrices Y, D, and X are constructed using static data, the proposed algorithm updates the dictionary matrices by conditioning dynamic observations (e.g., gas rate and bottomhole pressure (BHP)) to the ensemble.Figure 3 explains how to update the dictionary matrices during assimilations of the proposed algorithm.Let N para = N DCT and Y = Y DCT in the proposed algorithm.The weight matrix X ens is the state vector of the proposed algorithm.Only for the first assimilation, N ens realizations are randomly selected from X in order to compose the initial X ens .Thus, X ens is a N dict by N ens matrix while X is a N dict by N lib matrix.Let the superscript a denotes assimilation.Then, X a ens denotes the updated weight matrix after assimilation.Using D from Section 2.3 (Figure 3  ens .The quality is quantified in terms of ε qual , which indicates the discrepancy between observation and simulation results of an ensemble member: where N type is the number of data types. Let N qual be the number of qualified models selected from the ensemble.N qual models are selected among N ens models: Y a ens → Y a qual , as seen in Figure 3(d).Figure 4, for example, shows to select N qual = 20 realizations from N ens = 100 realizations regarding gas rate and BHP.It seems that the simulated responses of the selected models are closer to the reference results than the rest of the ensemble members.

Geofluids
After obtaining the qualified models, the library matrix Y DCT is updated as (Figure 3(f)) Note that ( 14) enhances both the quantity and quality of geologic libraries by conditioning Y DCT (which is conditioned to static observation data initially) to dynamic observation data with a relatively small number of reservoir simulation runs (N qual ≪ N lib ) for qualified ensemble members.
The updated Y DCT is decomposed to obtain newly updated D and X using OMP and K-SVD, as addressed in Section 2.3 (Figure 3(g)).With the new D and Y a DCT,ens in (11), X a ens is reupdated using OMP: Y a DCT,ens = DX a ens .This reupdated X a ens is used as m b in (2) as the state vector of the next assimilation.5 is a flowchart of the proposed ES-MDA algorithm that is compared with those of the other four ES-MDA algorithms investigated in this study.The general operating procedure for the algorithms is as follows: generation of an initial ensemble, reservoir simulation and selection of the qualified ensemble, Kalman-gain calculation and model update, facies designation, iteration of the above processes until the iteration number reaches the number of assimilations N a , and acquisition of qualified reservoir models.Differences between the algorithms are found in the types of state vectors and model update processes, as presented in Figures 5(b) and 5(c).

Framework of ES-MDA-DCT
Table 1 compares state vectors and sparse geologic dictionaries for the algorithms.ES-MDA updates gridblock facies indexes (originally assigned 0 and 1 for shale and sand, resp.).ES-MDA-DCT updates filtered DCT coefficients in the facies domain.ES-MDA-K-SVD tunes weights in facies domain with only one application of K-SVD before ES-MDA, while ES-MDA-DCT-K-SVD tunes weights in the DCT domain with only one application of K-SVD before ES-MDA.Finally, the proposed ES-MDA-DCT-i-K-SVD updates weights in the DCT domain with iterative K-SVD during ES-MDA.

Results and Discussions
The performance of the proposed algorithm (i.e., ES-MDA-DCT-i-K-SVD) is tested with applications to history matching of two channelized gas reservoirs.The algorithm performance is compared to those of four ES-MDA algorithms described in Table 1 to investigate coupling effects of dimensionality reduction and iterative sparse coding.Note that the developed algorithm updates dictionaries with dynamic data in each assimilation of ES-MDA.Neither ES-MDA nor ES-MDA-DCT adopts dictionaries.Dictionaries generated using static data are unchanged during data assimilation for ES-MDA-K-SVD and ES-MDA-DCT-K-SVD.
3.1.Field Description.Table 2 summarizes gas reservoir properties used in Case 1 and Case 2. For each case, the gas reservoir consists of sand and shale facies.Each facies has its relative permeability curves and absolute permeability value; permeabilities of sand and shale facies are 300 and 0.1, respectively.Reservoir boundaries are surrounded by numerical aquifers modeled by employing pore volume multipliers.
Figure 6 compares training images and reference models of the two cases.The size of the reservoir domain, well location, and well name are the same in both cases; eight vertical wells (P1, P4, P6, P7, P9, P12, P14, and P15) are drilled in the sand formation, and the other eight vertical wells (P2, P3, P5, P8, P10, P11, P13, and P16) are drilled in the shale formation.These 16 facies data at the well locations are regarded as hard data used for model generation.Both the reference model and N lib reservoir models are generated with hard data using SNESim that employs the training image shown in Figure 6(a) for Case 1.
The main differences between Case 1 and Case 2 are the shapes and widths of the channels.First, the reference model  3 describes operating conditions of 16 wells (from P1 to P16), the well coordinates of which are shown in Figure 6(a).Table 4 presents the number of parameters used in the five ES-MDA algorithms.Implementing DCT yields the data compression ratio N grid /N DCT = 12 1.All experiments were set up such that N a = 4 and α p = 4 according to (5).The proportion of the ensemble update in sparse geologic dictionaries N qual /N lib is 0.67%.7 shows the evolution of realizations in D achieved by invoking the proposed algorithm.We present five realizations having higher energy than the others, which are the first to the fifth column vectors in D, in each assimilation.Qualified realizations vary until the second assimilation is complete, implying a deficiency in the assimilation process due to a large dispersion of realizations.In the third assimilation, it seems that the first and fourth realizations conform to the reference model (Figure 6(a)).The first and second realizations in the fourth assimilation capture the S-shaped channel pattern of the reference model.This increase in the similarity among realizations implies an increase in the probability that sparse geologic dictionaries sufficiently represent the target gas reservoir.
Figure 8 compares history-matching results (day 1-day 3500) and prediction results (day 3501-day 7000) of the initial and updated ensembles against the reference results.Gas production rate and BHP are matching parameters of each ES-MDA algorithm.Water production rate is excluded from the matching dataset because it is a watch parameter.The three data types measured at wells P1, P4, P9, and P15 are shown because these wells installed in sand facies exhibit larger uncertainty than the other wells.In particular, water breakthrough does not occur at well P15 during the history-matching period.Thus, well P15 has the largest uncertainty in this case study.Solid gray, blue, dark blue, and red lines present well behaviors of the initial ensemble, the updated ensemble, the mean of the updated ensemble, and the reference model, respectively.The same initial ensemble is provided for each ES-MDA algorithm for a fair comparison.ES-MDA and ES-MDA-DCT expose a weakness in matching the reference results (Figures 8(a  11 Geofluids 16 wells as the breakthrough time is not detected during the matching period.It is reasonable to have moderate results at some matching targets despite the overall good matching performance [40,41]. The performances of the algorithms in terms of both history-matching (HM) and prediction (PD) periods are summarized in Table 5.The quality of the updated ensembles was assessed using the following: where μ ε is the normalized average of error ε (see (12)) and σ ε is the normalized standard deviation of ε.The superscripts upd and init refer to the updated ensemble member and the initial ensemble member, respectively.Smaller values of μ ε indicate a more accurate updated ensemble.Smaller values of σ ε indicate a smaller degree of reservoir uncertainty associated with the ensemble.Any value larger than 100% indicates deterioration of the ensemble quality, as revealed in the results obtained using ES-MDA and ES-MDA-DCT.In terms of gas rate, for example, the updated ensemble of ES-MDA-DCT amplifies the degree of dispersion compared to the initial ensemble.Furthermore, the behaviors of its ensemble mean do not follow those of the reference model.The proposed algorithm yields acceptable μ ε values and σ ε values that are smaller than those of the other ES-MDA algorithms.
Most observations are included within the bandwidths of the simulated profiles.The inclusion is also captured in the field cumulative production profiles of gas and water (Figure 9).Nonetheless, a remaining task is the construction of more robust dictionary matrices to improve the matching quality at every well target.Figure 10 compares ensemble mean maps and histograms of permeability obtained using the five ES-MDA algorithms.A task herein is to investigate whether the ES-MDA Table 5: Statistical parameters of history matching (HM) and prediction (PD) errors in terms of gas rate, water rate, and BHP (Case 1).μ ε and σ ε refer to the mean and standard deviation of ensemble error, respectively.12 Geofluids algorithms can capture the overall trend of the sand channel while preserving shale facies in the background.In the reference model, for example, the shortest path between wells P6 and P9 is filled with shale.Each ensemble member consists of sand facies with a permeability of 300 md and shale facies with a permeability of 0.1 md.However, the histogram of the initial ensemble mean (Figure 10 13 Geofluids that the unconformity is proportional to the frequencies of undesirable permeabilities between ln 0.1 (≈−2.3)md and ln 300 (≈5.7) md in the histograms.DCT is advantageous for preserving channel connectivity; however, an artificial sand channel connects wells P6 and P9.In addition, truncated DCT coefficients smooth out borders between the sand channel and background shale, as shown in Figure 10(c    conditioned to dynamic observations to the initially static library pool in each assimilation.Figure 13 compares production profiles obtained using the five ES-MDA algorithms.Each algorithm reuses the initial ensemble introduced in Case 1, as mentioned in Section 3.1.The proposed algorithm yields smaller standard deviations than the other algorithms and reduces reservoir uncertainty with reasonable matching accuracy (Table 7).The matching accuracy is comparable to those of the algorithms coupling sparse coding.Matches and prediction results for gas and water rate are decent.In terms of BHP, the proposed method shows the best conformance to observations except for well P9.The overestimated BHP is related to the underestimation of the water rate at well P9.Our inference is that this happens because of incorrect prior knowledge of the channel width.Field cumulative production profiles of gas and water are presented in Figure 14. Figure 15 compares averaged permeability distribution of the updated ensemble obtained using the five algorithms.The proposed algorithm outperforms the other algorithms by preserving the separation among the channels with NW-SE orientation and a distinct contrast at the borders between the shale background and sand channels.Facies index maps of individual ensemble members obtained after history matching also support the reliability of the proposed algorithm (Figure 16).

Conclusions
The hybridized ES-MDA algorithm coupled with iterative sparse coding in reduced parameter space was able to calibrate channelized reservoir models using and updating the repository of geologic features called sparse geologic dictionaries.The first library composed of thousands of reservoir models generated using SNESim was only conditioned to static data.Dimensionality reduction performed using DCT was effective in reducing the size of the library by converting gridblock facies into coefficients of discrete cosine functions.The weight matrix obtained by decomposing the library was imported to ES-MDA as a state vector.In each assimilation of ES-MDA, update of weights resulted in reservoir models that were well conditioned to both static and dynamic data.Adding the history-well-matched reservoir models as new members in the existing library was the critical factor for attaining improved matching accuracy and reduced model Table 7: Statistical parameters of history matching (HM) and prediction (PD) errors in terms of gas rate, water rate, and BHP (Case 2).μ ε and σ ε refer to mean and standard deviation of the ensemble error, respectively.dispersion because of the reflection of dynamic data in the updated dictionaries.The unprecedented consideration of dictionary update was the originality of this study and the contribution to researches about combinations of machine learning and ensemble-based data assimilation methods.
History-matching results of two channelized gas reservoirs (i.e., the S-shaped channel for Case 1 and three parallel channels for Case 2) indicated that the achievements arose from an iterative update of dictionaries of the proposed algorithm: the improved matching accuracy in both history and forecast in terms of well and total production, the reduced dispersion of production behaviors and permeability distribution, and the well-connected channel body of reservoir models with geological plausibility.ES-MDA with the  18 Geofluids dictionary update yielded higher matching accuracy values and lower dispersion values than ES-MDA incorporated with a fixed dictionary matrix.The increase in computational costs paid during the dictionary update was affordable compared to the assimilation algorithms not coupled with sparse coding.Improving the matching accuracy at some well-based levels remains as an outstanding task for the proposed technique despite the overall enhanced matching quality.We also anticipate that future works will adopt a more generalized sparse coding for diversifying the utility of the proposed framework in a variety of geoenvironments (e.g., naturally fractured reservoirs).
data measurement, processing, and interpretation.Stochastic characteristics of C D are reflected by z d ~N 0, I N d .z d is the random error matrix to observations, which is generated with a mean of zero and a standard deviation of I N d , where N d is the number of time steps in observations.

Figure 1 :
Figure 1: Example of discrete cosine transform (DCT) and inverse DCT (IDCT) applied to the reproduction of shale and sand facies of a channelized reservoir.
) and 2(b), each reservoir model consists of two facies: shale and sand with 0 and 1 indexes, respectively.Each column vector of Y corresponds to either a reservoir model or an encoded reservoir model because N para is determined depending on the type of assimilation algorithms in this study.N para = N grid (e.g., Y in Figure 2(b)) in the conventional ES-MDA, while N para = N DCT (e.g., Y DCT in Figure 2(c)) if DCT is applied to ES-MDA.In Figure 2(c), each column vector of Y DCT consists of DCT coefficients filtered from the corresponding column vector in Y.

Figure 2 :
Figure 2: Construction of sparse geologic dictionaries using DCT and K-SVD: (a) training image, (b) generation of initial channel models (Y) using SNESim, (c) transformation of Y into DCT coefficients, (d) construction of D and X from Y using K-SVD in DCT domain, and (e) reconstruction of Y from D and X.

Figure 2 (
Figure 2(d) illustrates how to decompose Y into D and X using K-SVD in the DCT coefficient domain.Note that column vectors in Y are sorted at random, while column vectors in D are sorted from left to right in descending order of energy to help us conveniently check the reservoir model quality.The term energy indicates the norm of a row vector in X corresponding to the realization in D. Greater norms in X indicate more essential principal components in D. The next task is how to decompose D and X from Y. Due to the large size of the three matrices, matrix decomposition is conducted numerically to minimize the discrepancy between Y and DX.In other words, min Y − Y′ F , 7 Figure 2(e) is an example to show that performing an IDCT yields Y DCT ′ while capturing geological features of the original Y DCT despite a diffusion in facies
(a)), Y a DCT,ens = DX a ens Updated ensemble in the first iteration 20 qualified members selected from the updated ensemble Mean of 20 qualified members

Figure 4 :
Figure 4: Gas rate and BHP profiles of 20 qualified realizations selected from the updated ensemble after the first assimilation.

Figure 5 :
Figure 5: Flow chart of five ES-MDA algorithms investigated in this study.The proposed algorithm corresponds to the number five in the boxes (b) and (c).

Figure 3 (
e) shows that facies of the N qual models are transformed into DCT coefficients: Y a qual → Y a DCT,qual , where Y a DCT,qual is the N DCT by N qual matrix composed of DCT coefficients of qualified ensemble members.

Figure 6 :
Figure 6: Training images and reference models adopted for history matching of (a) Case 1 and (b) Case 2.
) and 8(b)).Improvements in matching accuracy are accomplished in the simulation results of the updated ensembles obtained using ES-MDA coupled with sparse coding (Figures8(c)-8(e)).Executing sparse coding in the reduced parameter space has the result that most subfigures include the reference results within the bandwidth of the simulation results of the updated ensembles (Figures8(d) and 8(e)), thereby improving the matching accuracy.

Figure 7 :Figure 8 :
Figure 7: Evolution of realizations in the dictionary matrix D during multiple data assimilation (Case 1).

Figure 9 :
Figure 9: History matching results of Case 1: cumulative gas production and cumulative water production in the field.

Figure 10 :
Figure 10: Permeability distribution of the ensemble means and their histogram of permeabilities obtained after history matching of Case 1.
). ES-MDA-K-SVD, ES-MDA-DCT-K-SVD, and ES-MDA-DCTi-K-SVD reproduce channel patterns and the connectivity of the reference model (Figures 10(d)-10(f)), although Figure 10(e) contains undesirable sand facies between wells P6 and P9.A crisp contrast between shale and sand facies is observed in the outcome of the assimilation algorithms coupling DCT and K-SVD (Figures 10(e) and 10(f)).The developed algorithm reveals clearer shale facies in the path than the other algorithms (Figure 10(f)).

Figure 11 comparesFigure 11 :Table 6 :
Figure 11: Permeability distribution of five qualified reservoir models obtained after history matching of Case 1.

Figure 11 (
e) connects well P6 and well P9 with an artificial sand channel, but Figure11(f) does not.That is, conditioning dictionaries to dynamic data through the iterative sparse coding provides opportunities to calibrate the rock facies distribution between well P6 and well P9.The five realizations of the proposed algorithm keep the channel patterns similar to the reference, while maintaining shale facies on the shortest path between well P6 and well P9 (Figure11(f)).Most realizations have a broken sand channel near the upper left corner of the domain.This phenomenon would be an intrinsic limitation of geostatistical methods arising from a lack of observation data.Table6compares computational costs paid for the construction and update of dictionaries.The hardware specification utilized for comparison was Intel® Core™ i5-7200U @ 2.5 GHz with 8 GB RAM.The runtime of the reservoir simulation was excluded from the comparison.ES-MDA and ES-MDA-DCT cost nothing for dictionary construction.ES-MDA-K-SVD was a relatively costly algorithm.It took 218 minutes to obtain the original dictionary matrices D and X via Y despite the matrix construction one time.DCT helped save computational costs required for sparse coding.ES-MDA-DCT-K-SVD was the cheapest algorithm (taking 5.7 minutes) thanks to the reduced number of parameters due to DCT.Updating sparse geologic dictionaries in each assimilation resulted in an additional 15.7 minutes in the subsequent four assimilations.Nevertheless, it seems this amount of extra cost is within a reasonable range.ES-MDA-DCT-K-SVD was 38.2 times as fast as ES-MDA-K-SVD, and ES-MDA-DCT-i-K-SVD was 10.2 times faster than ES-MDA-K-SVD.

Figure 15 :
Figure 15: Permeability distribution of the ensemble means and their histogram of permeabilities obtained after history matching of Case 2.
Covariance matrix of state vectors C D : Covariance matrix of observed measurement error C dd : Autocovariance matrix of simulation data d C md : Cross-covariance matrix of state vector m and simulation data d d: Simulation data d obs : Observation data d unc : Perturbed observation data d r : Column vector of the dictionary matrix D (i.e., realization) d: Mean of simulation data D: Dictionary matrix E j : Discrepancy between the library matrix Y and the reconstructed library matrix Y ′ except the jth d r j x j f : Reservoir simulator J: Objective function m: State vector of a reservoir model m: Mean of state vectors m b : State vector of a reservoir model before update I N d : N d by N d identity matrix

Figure 16 :
Figure 16: Permeability distribution of five qualified reservoir models obtained after history matching of Case 2.

Table 1 :
State vectors and update processes of sparse geologic dictionaries for five ES-MDA algorithms investigated in this study.

Table 2 :
Experimental setting of reservoir parameters used in Case 1 and Case 2.
tion.Secondly, the channel width of the training image of Case 2 is 1.2 times thicker than that of Case 1, making that the channel width of the reference model of Case 2 is 1.2 times thicker as well.Despite the differences, this study intentionally reuses the initial ensemble designed for Case 1 as the initial ensemble of Case 2. The manipulation of an initial

Table 3 :
Experimental data for well parameters used in Case 1 and Case 2.

Table 4 :
Number of parameters used for construction and update of sparse geologic dictionaries.