Next Article in Journal
Abnormal Data Cleaning Method for Wind Turbines Based on Constrained Curve Fitting
Next Article in Special Issue
Development of Virtual Flow-Meter Concept Techniques for Ground Infrastructure Management
Previous Article in Journal
Solid-Rotor Induction Motor Modeling Based on Circuit Model Utilizing Fractional-Order Derivatives
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Joint History Matching of Multiple Types of Field Data in a 3D Field-Scale Case Study †

by
William Chalub Cruz
1,*,
Xiaodong Luo
2 and
Kurt Rachares Petvipusit
3
1
University of Stavanger, 4021 Stavanger, Norway
2
Norwegian Research Centre (NORCE), 5008 Bergen, Norway
3
Equinor, 5254 Bergen, Norway
*
Author to whom correspondence should be addressed.
This paper is modified from the material presented in SPE Norway Subsurface Conference, Bergen, Norway, 27 April 2022.
Energies 2022, 15(17), 6372; https://doi.org/10.3390/en15176372
Submission received: 2 August 2022 / Revised: 27 August 2022 / Accepted: 28 August 2022 / Published: 31 August 2022
(This article belongs to the Special Issue Advances in Reservoir Simulation)

Abstract

:
This work presents an ensemble-based workflow to simultaneously assimilate multiple types of field data in a proper and consistent manner. The aim of using multiple field datasets is to improve the reliability of estimated reservoir models and avoid the underestimation of uncertainties. The proposed framework is based on an integrated history matching workflow, in which reservoir models are conditioned simultaneously on production, tracer and 4D seismic data with the help of three advanced techniques: adaptive localization (for better uncertainty quantification), weight adjustment (for balancing the influence of different types of field data), and sparse data representation (for handling big datasets). The integrated workflow is successfully implemented and tested in a 3D benchmark case with a set of comparison studies (with and without tracer data). The findings of this study indicate that joint history matching using production, tracer and 4D seismic data results in better estimated reservoir models and improved forecast performance. Moreover, the integrated workflow is flexible, and can be extended to incorporate more types of field data for further performance improvement. As such, the findings of this study can help to achieve a better understanding of the impacts of multiple datasets on history matching performance, and the proposed integrated workflow could serve as a useful tool for real field case studies in general.

1. Introduction

In the petroleum industry, reservoir simulation plays an important role for the decision-making process of oil field developments, where a common practice is to apply a closed-loop workflow [1] to develop and manage reservoir performance under different scenarios and conditions, as shown in Figure 1. This type of simulation is commonly used to predict and/or optimize future oil, gas, and water productions in a field (referred to as system). However, uncertainties and simplifications are involved in the construction of reservoir models (system models), due to the lack of sufficient information in, e.g., geology, seismic data, well logs and tests. Consequently, reservoir parameters (e.g., permeability, and porosity) are not accurately obtained, which makes reservoir models uncertain for production forecasts, and the corresponding decision-making process unreliable for reservoir management. In order to improve the reliability of production forecasts, one needs to condition reservoir models on sensored field datasets, through a process known as history matching in the literature. In this way, the estimated parameters of reservoir models can generate simulated outputs that are as close as possible to the field datasets.
In a standard history matching process, the field data used to estimate uncertain reservoir model parameters are well production data, which include bottom hole pressure, and oil, gas, and water rates. These kinds of measurements are frequent in time and relatively cheap in terms of acquisition cost, but their usefulness is often somewhat insufficient for the purpose of reducing uncertainties in production forecasts, due to their limited spatial coverage of the reservoir [2,3]. In addition to production data, other types of field data, e.g., inter-well tracer and 4D seismic data, can be used as complementary sources of information to further reduce the uncertainties-hence improve the reliability-of reservoir models, which constitutes the focus of the current study.
Inter-well tracer test (IWTT) has been proven as an efficient technology to obtain information of fluid dynamics, well-to-well communication, and heterogeneities such as fractures and flow barriers (for a review of tracer studies see [4]). In the tracer technology, inert compounds (e.g., radioactive, chemical, or natural tracers) are used to label water or gas from specific wells, and trace fluid movements as the tracer flow moves through the reservoir together with the injection fluid. After the first breakthrough in a producer, IWTT provides a reliable and definite information of a well-to-well communication.
On the other hand, 4D seismic data contain much more information in space, especially about pressure and fluid saturation changes. Therefore, seismic data provide another useful source of information for understanding reservoir behavior, and identifying zones of remaining oils [5]. When used in combination with production and tracer data, they could further improve the qualities of reservoir models through history matching, and thus generate more reliable production forecasts for better decision making.
Despite the appealing features of tracer and 4D seismic data, they are still often qualitatively employed for reservoir monitoring and management in the petroleum industry, but are relatively less used in a quantitative way, due to various reasons arising from both theoretical understandings and practical applications to jointly match both types of field data in a coherent way. Tracer has been incorporated into history matching by a few authors [3,5,6,7,8,9,10]. Meanwhile, recently, several authors have demonstrated the benefits of incorporating seismic data into a history matching process [2,11,12,13,14,15]. Nevertheless, methodologies for jointly history matching production, tracer, and 4D seismic data have not been previously studied in the literature yet. Thus, there is a potential to include both tracer and 4D seismic data into a history matching workflow to further reduce the uncertainties and improve the reliability of reservoir models.
In the literature, there are different schools of approaches to tackling history matching problems. For instance, history matching can be formulated as a conventional minimization problem, in which an iterative algorithm is used, in combination with certain geoestatistical modelling methods [16,17], to find the best matched model(s). Among the iterative history matching algorithms, ensemble-based methods have received immense attention in the petroleum industry, for their convenience in implementation and computational efficiency to handle big models and big data sets when compared to other conventional methods [18]. For instance, the ensemble Kalman filter (EnKF) introduced by [19] was initially extensively used for history matching problems [20,21]. As a sequential algorithm, however, the EnKF requires frequent simulator restarts, which may cause significant additional simulation time and other practical challenges (e.g., more prone to failed reservoir simulations) [22]. To avoid these noticed issues, ensemble smoother (ES) [23] and its iterative versions [24,25,26,27,28] can be adopted instead. In comparison to the EnKF, ES and iterative ES (IES) do not need simulator restarts and have less model variables to update during history matching, and are practically easier to implement than the EnKF for history matching problems. In addition, various numerical results (e.g., [25,26]) indicate that the IES tends to perform better than the EnKF and the ES. As a result, currently, the IES appears to be among the state-of-the-art approaches to history matching problems, and is thus adopted in the current work.
The recent developments of ensemble-based history matching algorithms make the quantitative use of multiple field data sets easier and faster. Despite this technological advancement, to the best of our knowledge, there are no existing studies demonstrating the benefits of simultaneously assimilating production, tracer and 4D seismic data into reservoir models through an ensemble-based history matching algorithm. As a main focus (and contribution) of this study, our primary objective is to demonstrate the benefits of assimilating tracer and 4D seismic data (in addition to production) in an integrated, ensemble-based history matching workflow, through a 3D field-scale case study.
This work is organized as follows: First, we provide details of the integrated ensemble-based history matching workflow, which includes the formulation of an IES algorithm and a correlation-based adaptive localization scheme. Second, we apply the history matching workflow to a 3D field case, the Brugge benchmark [29], and examine the performance of the proposed workflow. Finally, we conclude the work with some technical discussions and our future research plan.

2. Ensemble-Based History Matching Workflow

Figure 2 illustrates the flowchart of the integrated history matching workflow using multiple types of field data. To formulate the history matching problem, we assume there is a forward simulator g (e.g., a numerical reservoir and/or seismic simulator) which outputs a N d -dimensional vector containing the simulated data d s i m , given a N m -dimensional vector of reservoir model m as the input:
d s i m = g ( m ) .
Note that, in Step 1 of Figure 2, the forward simulator generates an ensemble of simulated production and tracer data, with an ensemble of input reservoir models, while the forward seismic simulator produces an ensemble of simulated amplitude-versus-angle (AVA) data, with the associated petro-elastic model (PEM) and AVA at Steps 2 and 3, respectively, based on the simulated pressure and saturation profiles from the forward reservoir simulator, and the simulated velocities and density profiles from the PEM, respectively. The formulation of the forward seismic simulator is presented in Appendix A.
Furthermore, we have the observed (field) data d o which are obtained through the following noisy observation system:
d o = g ( m true ) + δ ,
where m true stands for the ground-truth model (unknown in a real reservoir), g ( m true ) for the true output of the forward model, and δ for a N d -dimensional vector of contamination noise that may be present in the course of field data acquisition, which are assumed to follow a multivariate Gaussian distribution with zero mean and covariance (referred to as measurement error-covariance matrix) represented by C d , i.e., δ N ( 0 , C d ) .
With the above observation system, in Step 4 of Figure 2, a history matching algorithm is then adopted to find one or multiple reservoir models m whose simulated output g ( m ) matches the observed data d o reasonably well. In practice, when the size of the observed data is much smaller than of the reservoir model parameters (e.g., N d N m ), history matching is an under-determined and high-dimensional inverse problem, which makes the history matching problem challenging, in the sense that in principle there could exist an infinite number of reservoir models matching the observed data equally well (non-uniqueness), due to the large degree of freedom (DOF) from the model side [30].
To mitigate this problem, one can either introduce a certain regularization term in the form of a least-square problem, or increase the number of field data in history matching [3]. Here we consider both aforementioned ideas, by including both a regularization term into a relevant cost function (cf. Equation (3)), and multiple types of field data into history matching.
Practical challenges may arise with multiple types of field data in history matching. For instance, different from production and tracer data, 4D seismic data are less frequent in time, but with a much larger size. Thus, when assimilating large datasets, computational issues regarding storage memory and CPU time may emerge. Meanwhile, large datasets could have a dominant influence on model updates in history matching. Moreover, with large field datasets, ensemble-based history matching methods are often prone to ensemble collapse (meaning that an ensemble of reservoir models also collapses into a single one, with few varieties among reservoir models), and thus tend to under-estimate model uncertainties [15]. Our approaches to handling the aforementioned issues include conducting a sparse representation of seismic data for dimension reduction, adjusting the relative weights among different types of field data to balance their influence on history matching, and conducting correlation-based localization to mitigate the issue of under-estimated uncertainties, which will be elaborated in Section 3 later.

2.1. Iterative Ensemble Smoother Based on a Regularized Levenburg–Marquardt Algorithm

As history matching problems are typically non-linear, a certain iterative method is needed to deal with the non-linearity. In the current work, we adopt one such method, called iterative ensemble smoother based on the regularized Levenburg–Marquardt algorithm (IES-RLM) [26].
In IES-RLM, one aims to find an ensemble M i + 1 { m j i + 1 } j = 1 N e of reservoir models that approximately minimizes a nonlinear-least-squares (NLS) cost function, as follows:
arg min { m j i + 1 } j = 1 N e 1 N e j = 1 N e L j i + 1 m j i + 1 , L j i + 1 m j i + 1 d j o g m j i + 1 T C d 1 d j o g m j i + 1 + γ i m j i + 1 m j i T C m i 1 m j i + 1 m j i ,
where i stands for the index of iteration step, j for the index of ensemble member, N e for the ensemble size, and C m i for the sample model error-covariance matrix with respect to the prior ensemble M i { m j i } j = 1 N e , in the form of C m i = S m i ( S m i ) T , with the squared-root matrix S m i defined in Equation (5).
The NLS-type cost function in Equation (3) contains two terms: the first one (called data mismatch term) calculates the difference between observed and simulated data, and the second one (called regularization term) imposes a constraint that aims to prevent a large deviation of the updated model m j i + 1 from its predecessor m j i . Here, the regularization form is known as Tikhonov regularization [30], but one can also use other types of regularization [27]. Meanwhile, γ i > 0 in Equation (3) is a coefficient influencing the relative weight between the data mismatch and regularization terms, whose value changes over the iteration steps [26]. In this optimization problem, the goal is to minimize the average of an ensemble of cost functions at each iteration step, until a certain stopping criterion is reached. The stopping criteria used in this work consist of the following three ones: (1) the maximum number of iteration is reached; (2) the relative change of data mismatch values in-between two consecutive iteration steps is less than 0.01%; (3) the data mismatch value at a given iteration step is lower than the number of observed data points (Nd) times a factor (1 in this study).
By minimizing the cost function in Equation 3, one can obtain the following approximate model update formula:
m j i + 1 = m j i + S m i S ˜ d i T S ˜ d i S ˜ d i T + γ i I 1 Δ d ˜ j i ,
where the square root matrix S m i in the model space is defined as
S m i = 1 N e 1 m 1 i m ¯ i , , m N e i m ¯ i ;
m ¯ i = 1 N e j = 1 N e m j i .
Similar to Equation (5), one can define a square root matrix S d i in the data space as
S d i = 1 N e 1 g m 1 i g m ¯ i , , g m N e i g m ¯ i .
In many practical history matching problems, the observations may contain different types of field data in different orders of magnitudes. To mitigate potential issues related to these imbalanced magnitudes, one can introduce a normalization procedure to quantities in the observation space, so that S d i and d j o g m j i are normalized by a square root C d 1 / 2 of C d ( C d = C d 1 / 2 ( C d 1 / 2 ) T ) to arrive at the following formulations:
S ˜ d i = C d 1 / 2 S d i ,
Δ d ˜ j i = C d 1 / 2 d j o g m j i .
In a practical implementation of the model update formula Equation (4), numerical stability could be an issue when inverting the matrix ( S ˜ d i S ˜ d i T + γ i I ) (especially when using 4D seismic data), since the data size N d is typically much larger than the ensemble size N e . Therefore, a common procedure introduced in the literature to obtain a numerically more stable algorithm is to carry out an inversion for a certain matrix with a lower dimension (less than the ensemble size N e ) by applying a Truncated Singular Value Decomposition (TVSD) to the matrix S ˜ d i [24,31]. Suppose that, through the TSVD, one obtains
S ˜ d i U ^ i Σ ^ i ( V ^ i ) T ,
where U ^ i [ u 1 i , , u sv i ] and V ^ i [ v 1 i , , v sv i ] are unitary matrices consisting of the kept left and right eigenvectors of S ˜ d , respectively, and Σ ^ R sv × sv is a rectangular diagonal matrix containing in its diagonal a set of retained leading singular values. The integer sv is chosen as suggested by [31] to keep the leading singular values that add up to at least 99% of the total sum of squared singular values.
Inserting Equation (10) into Equation (4), one obtains a modified update formula
m j i + 1 = m j i + S m i V ^ i ( Σ ^ i ) T Σ ^ i ( Σ ^ i ) T + γ i I 1 ( U ^ i ) T Δ d ˜ j i .
In the current work, the regularization parameter γ i in Equation (11) is chosen as
γ i = η i × T r ( Σ ^ i ) T Σ ^ i / N sv ,
where η i is a positive number that starts with a value of 1 and varies as a function of the iteration step, following the rule in [26]; and the operator T r calculates the trace of a matrix.

2.2. Correlation-Based Automatic and Adaptive Localization

In practical applications of ensemble-based methods, a limited number of reservoir models is typically adopted to reduce computational cost. Under this setting, certain numerical issues may arise when the number of reservoir models N e is significantly smaller than the sizes of reservoir model N m and field data N d . For instance, in the context of ensemble-based history matching, the limited number of reservoir models can produce spurious correlations between uncorrelated reservoir model parameters and observation data, leading to unsatisfactory history matching performance due to the excessive reduction of ensemble variability (e.g., ensemble collapse). To mitigate the negative effects of the limited number of reservoir models, it is common to adopt an auxiliary technique, called (Kalman gain) localization [32,33,34,35] for ensemble history matching algorithms.
For an easier explanation of the main idea behind localization, one can rearrange Equation (11) as
m j i + 1 = m j i + K i Δ d ˜ j i ,
with the Kalman gain matrix K i R N m × N d defined as
K i = S m i V ^ i ( Σ ^ i ) T Σ ^ i Σ ^ i T + γ i I 1 ( U ^ i ) T .
To conduct localization in Equation (13), one can replace the Kalman gain K i by a Schur product (element-wise product) between K and a localization matrix C , as in
m j i + 1 = m j i + C K i Δ d ˜ j i .
Here, the localization matrix C is implemented as a tapering matrix, introduced to assign different weights (tapering coefficients) for different combinations of reservoir model parameters and field data points. Note that, in its conventional form, the Kalman gain matrix needs to be computed at each iteration step, which could be a very high-dimensional matrix, especially with big observation data (e.g., 4D seismic data) and big reservoir models. To deal with the high dimensionality of K i , one can choose to sparsely represent big observation data in another domain, so as to obtain a much smaller number of data points used in history matching, as will be explained later.
Equivalently, Equation (15) can be expressed as:
m j , k i + 1 = m j , k i + s = 1 N d ( c ks k ks i ) Δ d ˜ j , s i = m j , k i + s = 1 N d k ks i ( c ks Δ d ˜ j , s i ) .
where c ks [ 0 , 1 ] and k ks i stand for the elements of C and K i on the k-th row and the s-th column, respectively, Δ d ˜ j , s i is the s-th element of Δ d ˜ j i , and m j , k i + 1 , and m j , k i the k-th element of the vectors m j i + 1 , and m j i , respectively.
There are different ways to compute the tapering coefficients, c ks , in the literature. Among them, a common approach appears to be distance-based localization [32,33,36], in which the observations and reservoir model variables are assumed to have physical locations, so that one can compute the distance between the physical locations of a model variable and a field data point. This approach can work well in many situations. However, there are some notable issues. For instance, the observation data need to have associated physical locations, the tapering function is not adaptive to reservoir heterogeneities, and the difficulty to “localize” non-local observations when there is a long range of correlations between model variables and field data.
To deal with the aforementioned issues, ref. [35] proposed a correlation-based adaptive localization scheme. In this approach, the authors computed the tapering values c ks dependent on the sample correlation ρ k s between the k-th element of the initial ensemble m j , k 0 and the corresponding initial ensemble of the s-th innovation element Δ d ˜ j , s 0 ( j = 1 , 2 , , N e ). Then, they defined a threshold value θ s dependent on the noise level (standard deviation, STD) of the sampling errors. To compute the noise level, one can assume that the members of the initial ensemble m j 0 ( j = 1 , 2 , , N e ) are independent and identically distributed (i.i.d). Due to the independence assumption, one can obtain a new ensemble m ^ j 0 by randomly shuffling the indices j of m j 0 . After that, one can compute the sampling errors ρ ^ s , which represents the correlation field between the new ensemble m ^ j 0 ( j = 1 , 2 , , N e ) and the ensemble of the s-th innovation-data point Δ d ˜ j , s i ( j = 1 , 2 , , N e ) . Under the i.i.d assumption, m ^ j 0 and m j 0 are independent. As a result, one can take the correlation field ρ ^ s as a realization of sampling errors of the sample correlations between all model parameters and the s-th innovation-data point. Then, for each type of petro-physical parameter and each field data point, one can estimate the noise level σ ^ s with respect to the sampling errors ρ ^ s as follows [35]:
σ ^ s = median ( abs ( ρ ^ s ) ) 0.6745 .
After that, one can further compute the threshold value θ s by
θ s = σ ^ s 2 ln ( # ρ ^ s ) ,
where # ρ ^ s is the number of elements in ρ ^ s . Note that one should perform this procedure for each type of petro-physical parameter and each field data point.
Finally, one can generate a smooth tapering function as in [35]
z = 1 abs ( ρ k s ) 1 θ s ,
where the k-th element ρ k s of ρ s is the sample correlation between a k-th model variable and the s-th field data point. Then, the variable z is used in the Gaspari–Cohn function (Gaspari and Cohn, 1999)
c ks = f G C ( z ) = 1 4 z 5 + 1 2 z 3 + 5 8 z 3 5 3 z 2 + 1 , if 0 z 1 1 12 z 5 1 2 z 4 + 5 8 z 3 + 5 3 z 2 5 z + 4 2 3 z 1 , if 1 < z 2 0 , if z > 2 ,
to finally obtain the tapering coefficient c ks .

3. Application to the Brugge Field Case Study

The Brugge field model is a benchmark case based on a North sea field and developed by [29], and it has the characteristics and complexities of reservoir models used in real field case studies. The Brugge field consists of two zones, one with higher permeability and porosity, and the other with lower permeability and porosity. The zone with high permeability and porosity is located on layers 1–2 and layers 6–8, while the other with low permeability and porosity is located on layers 3–5 and 9.
The dimensions of the reservoir model are 139 × 48 × 9, summing up to 60,048 gridblock, in which 44,550 are active. Figure 3 shows the location of the field wells, including 30 wells (20 inner producers and 10 injectors). The producers are located more centrally in the red area and labelled as BR-P-1-BR-P-20, while the injectors are located on the border and labelled as BR-I-1-BR-I-10, with water being the only injected fluid.
The initial ensemble of this benchmark contains 104 realizations of reservoir models, from which we take the first realization (#1) as our reference (true) case to generate the observed data. As for reservoir model variables, we consider porosity (PORO) and permeability along x-, y-, and z- directions (denoted by PERMX, PERMY, and PERMZ), summing up to 178,200 ( 4   × 44,550) uncertain petro-physical parameters to estimate in history matching.
For illustration, Figure 4 and Figure 5 show the distributions of the porosity and the log permeability (along the x-direction), respectively, on each layer from realization number two ( # 2 ) of the initial ensemble. As reported in [29], the permeability and porosity are linearly correlated, and the spatial distribution of the permeability is anisotropic, meaning that permeability distributions along different directions may not be the same.
In history matching, we condition reservoir models on three types of field data, namely, production, tracer, and 4D seismic data. The historical production data cover a period of 3647.5 days, with production forecasts from 3647.5 to 9869.5 days, where reservoir simulations are conducted using a black oil simulator (ECLIPSE©). The production data include well bottom hole pressure (WBHP), well oil and water rates (WOPR, and WWPR), summing up to 1400 data points. In addition, an injection of a water passive tracer pulse of 1.0 lb per STB is performed from day 749 to day 1812 in the injector well BR-I-6 (WTPCW06), summing up to 400 data points. To mimic the presence of measurement noise in real production data, here we also introduce certain zero-mean Gaussian white noise to the reference production and tracer data, for which the STD of the noise is set as 10% of the magnitude of each production or tracer data point, except for WBHP (for which the noise STD is set as 50 psi for each data point instead).
The 4D seismic data are obtained from three surveys: base (day 1), monitor #1 (day 991), and monitor #2 (day 2999). From each seismic survey, the seismic attributes (AVA data) are obtained using two different offset angles, near (10 ) and mid (20 ). The AVA data are in a seismic domain and do not possess physical locations in the reservoir coordinate system, meaning that the AVA data are in a different dimension compared to the dimension of the reservoir model. In this work, the dimension of the AVA data is 139 × 48 × 176 for each seismic survey, summing up to 7,045,632 data points in total [15]. We also introduce certain zero-mean Gaussian white noise to the AVA data, with the STD of the noise set to 30% of the magnitude as in [15].
In order to deal with the issue of big seismic data in history matching, and thereby the required memory, we choose to sparsely represent the AVA data by first applying 3D discrete wavelet transform (DWT) to the seismic data, and then (through a thresholding operation) selecting a small set of the resulted leading wavelet coefficients as the observations in history matching. For brevity, we skip the details of the DWT-based sparse data representation procedure. Readers are referred to [15,37] for more information.
From the algorithmic perspective, in the presence of the DWT-based sparse data representation procedure, the effective observation system with respect to 4D seismic data becomes
T W ( d o ) = T W g s ( m true ) + W ( δ ) ,
where W and T denote wavelet transformation and thresholding operators, respectively, and g s corresponds to the seismic forward simulator, as discussed in Appendix A. By using the aforementioned formulation, we are able to represent the original 4D seismic (7,045,632 data points) by a much smaller set of 1896 wavelet coefficients, as elaborated in [15].

Experiment Settings

In our previous work [3], we have shown that adopting both production and tracer data in the Brugge benchmark helps improve the performance of history matching, in comparison to the choice of using production data only. The current work can be considered as a follow-up study, in which we aim to further examine the impacts of 4D seismic data on history matching, and show the complexity of the joint history matching problem in the presence of multiple types of field data. To this end, in the current work we conduct two experiments complementary to those in [3]:
  • Case 1: using production and 4D seismic data;
  • Case 2: using production, tracer, and 4D seismic data.
In both cases, we apply the DWT-based sparse data representation procedure to the seismic attributes obtained at three surveys for improved computational efficiency.
We use IES-RLM as the history matching algorithm, and equip it with the correlation-based adaptive localization scheme described in Section 2.2. We start the IES-RLM algorithm with η 0 = 1 , with the reduction and increment factor being, 0.9 and 2, respectively. More specifically, if the average data mismatch is reduced at the current iteration step, then the η value at the next iteration step is set to 0.9 times that at the current step; otherwise, the next η value is set to 2 times the value of the current iteration step, and an inner loop is then activated to search for lower data mismatch.
To simultaneously history match multiple types of field data, we introduce a normalization procedure to adjust the relative weights assigned to different types of field data, as in [2]. Following the formulation in [3], we substitute the measurement error-covariance matrix C d in Equation (3) by the Schur product diag w C d , where w is a vector to be specified later, and the operator diag w stands for a diagonal matrix whose diagonal elements are taken from the vector w .
For production and tracer data, the elements of w are determined using the following rule:
w k = m a x ( type ( d k o ) ) × p ( type ( d k o ) ) / σ k 2 , for k = 1 , 2 , , N d ,
where the operator m a x ( type ( d k o ) ) takes the maximum value of a given type of field data, and type ( d k o ) stands for the type (WBHP, WOPR, WWPR or WTPCW06) of the k-th field-data point. Note that the value m a x ( type ( d k o ) ) is normalized by the measurement-noise STD σ k associated with the data point d k o , aiming to mitigate potential problems caused by the different orders of magnitudes of the field data. The notation p represents a percentage value, which is set to 5 % if d k o is a type of production data (WBHP, WOPR, or WWPR), and to 20 % instead if d k o is a tracer data point (WTPCW06).
For seismic data, the elements of w are determined using the following rule:
w k = m e d i a n ( a b s ( W ^ ( d k o ) ) ) 0.6745 / σ k 2 ,
where W ^ ( d k o ) is the wavelet coefficient belonging to the finest sub-band of W ( d k o ) .
As aforementioned, 4D seismic data usually contain a large number of data points, which can dominate the updates during history matching. To overcome this noticed issue, a certain scaling factor s is introduced to adjust the relative weights assigned to production, tracer and 4D seismic data during history matching. The scaling factor is determined using data mismatch of the initial ensemble, which is computed using the following formula:
HM j i = ( Δ d ˜ j i ) T Δ d ˜ j i .
By using Equation (24), one can compute the scaling factor as follows:
s = H M ¯ 1 0 / H M ¯ 2 0 ,
where the overline denotes the mean value over the ensemble members, H M ¯ 1 0 corresponds to the mean data mismatch value of 4D seismic data, and H M ¯ 2 0 to that of production data (Case 1) or that of both production and tracer data (Case 2). Note that in Case 2, we group production and tracer data together to compute the mean of data mismatch. To construct the scaled observational data, one needs to multiply the observation data from group 2 by the scalar s, so that Δ d ˜ j , 2 i is replaced by [ Δ d ˜ j , 2 i × s ] .
To assess uncertainty reduction for each type of petro-physical parameter (PERMX, PERMY, PERMZ, and PORO) in the experiments, we use the Sum of Normalized Variance (SNV) [38], as in
SNV = k = 1 N m var ( m k f ) var ( m k 0 ) ,
where var ( m k 0 ) and var ( m k f ) denote the variances of a particular type of petro-physical parameter distributed on the k-th active reservoir gridblock of the initial and final ensembles, respectively.
Finally, to assess and compare history matching performance, in terms of the discrepancy between an estimated reservoir model m and the true one m true , we use the Root Mean Squared Error (RMSE), as in
RMSE = m m true 2 N m .

4. Results and Discussions

This section focuses on illustrating the performance of the history matching of the two aforementioned experiments, in terms of data mismatch values during history matching (which is calculated through Equation (24)) and forecast periods, and also uncertainty quantification that is represented by mean, standard deviation, and RMSE of the final ensemble.
Table 1 summarizes the values of data mismatch ( m e a n ± S T D ) and RMSE ( m e a n ± S T D ) with respect to the initial ensemble and the final ensembles of the two experiments’ settings (Cases 1 and 2). In comparison with the initial ensemble, both experiments exhibit better history matching performance, in terms of both lower data mismatch and RMSE values. Comparing the results between the two different experiments shows that we achieve slightly lower data mismatch for production data (−0.2389%), but at the cost of a slightly higher data mismatch for 4D seismic data (+0.9862%), when jointly history matching production, tracer, and 4D seismic data (Case 2). In terms of RMSE, Case 2 contains lower average values for PERMY and PERMZ (−1.0565% and −3.3100%, respectively), but slightly higher ones for PERMX and PORO (+1.6433% and +0.7905%, respectively). In addition, the average RMSE values for all parameters together in Case 2 is slightly lower in comparison with Case 1, which indicates that the overall history matching performance in terms of RMSE tends to improve by adding more types of field data, but the complexity of the history matching process tends to increase. Note that in Case 2, only inter-well tracer from injector BR-I-6 was included in the history matching workflow, which provides additional information regarding well-to-well connectivity for the producers (BR-P-3, BR-P-4, BR-P-14, BR-P-15, BR-P-16, and BR-P-17), where tracer is detected. However, one cannot expect a substantial reduction in the overall data mismatch and RMSE values due to the use of the tracer data, as the inter-well tracer is usually detected in a limited number of production wells in the reservoir and tracer from a single injector well may not be sufficient to provide information to improve the connectivity of the entire reservoir. We discuss more about reservoir connectivity in the forecast analysis.
For illustration, Figure 6 reports the ensemble of data mismatch values obtained at each iteration step, in the form of box plots in Cases 1 and 2. In both cases, data mismatch values tend to decrease over the iteration steps for all types of field data. In addition, one can observe that main reductions of data mismatch values take place within the first few iteration steps for production data, and afterwards the changes of data mismatch values are less substantial. These main reductions in data mismatch are related to the changes of well bottom hole pressure (WBHP), which is strongly reduced in the first few iterations steps. Furthermore, the IES-RLM algorithm is close to convergence after 10 iterations steps, which is the maximum number of iterations used in the case studies.
Similarly, Figure 7 shows box plots of RMSE values in the two case studies, for PERMX, PERMY, PERMZ and PORO, respectively. As one can see, the average RMSE values also tend to decrease over the iterations steps. Consistent with Table 1, the results show that jointly history matching production, tracer and 4D seismic data (Case 2) tends to result in lower mean RMSE values for PERMY and PERMZ, but slightly higher ones for PERMX and PORO, in comparison with the choice of history matching both production and 4D seismic data (Case 1).
Figure 8 presents the profiles of WOPR and WWPR from well BR-P-15 and BR-P-16, with respect to the initial and final ensembles in the two case studies. In comparison to the initial ensemble, both cases indicate significant improvements in terms of reduced data mismatch and uncertainty (in terms of ensemble variability). For instance, compared to the choice of history-matching production and 4D seismic data (Case 1), one can see that the inclusion of tracer in Case 2 tends to slightly reduce the variability of the final ensemble of WWPR in BR-P-16, which means that the final ensemble follows the observed data better.
Similar to Figure 8, Figure 9 shows the profiles of tracer data injected in well BR-I-6 and recorded at wells BR-P-15 and BR-P-16, with respect to the initial and final ensembles in the two case studies. Here, the history matching performance is reasonably good in terms of data mismatch, e.g., for both production wells, history matching helps reduce data mismatch of the final ensembles, which is consistent with Table 1.
In addition, Table 2 provides an assessment of parameter uncertainty reduction (in terms of SNV) in the two aforementioned experiments. As one can see, in both case studies, the final ensembles maintain substantial variability, meaning that the localization scheme is useful for preserving ensemble variability. Consistent with Table 1, Case 2 contains slightly lower SNV values for PERMY and PERMZ, while slightly higher ones for PORO and PERMX, in comparison with Case 1. As tracer data are more correlated to permeability, one may expect more uncertainty reduction for permeability than porosity. This is partially reflected by the results in Table 2, where there is a relatively stronger uncertainty reduction for PERMY and PERMZ (−0.51% and −3.09%, respectively).
Another important aspect is to inspect the spatial distributions and standard deviations of the petro-physical parameters after history matching. Figure 10 and Figure 11 report mean and standard deviation maps of the petro-physical parameters on layer 2, with respect to the initial and final ensembles in Cases 1 and 2. After history matching, one can see that there are minor changes among different experiment results. Relative to the mean maps, our history matching workflow is able to capture the main aspects of the initial ensemble in both experiments. Meanwhile, the standard deviation maps of the two experiments indicate substantial variability after history matching, which means that the ensemble collapse is avoided.
Finally, we forecast the production data beyond the history matching period using the final ensembles obtained in Cases 1 and 2. For comparison, we calculate the forecast data mismatch for WOPR and WWPR, which is the difference between simulated data of the reference model and an estimated ensemble of reservoir models at the forecast period, normalized by the STD of the WOPR and WWPR.
As aforementioned, we observed better history matching performance in terms of production data mismatch and RMSE values for PERMY and PERMZ parameters, when using production, tracer, and 4D seismic data (Case 2) during the history matching period. In terms of forecast mismatch, the comparison becomes a bit more complicated. As one can see in Figure 12, in some wells (e.g., BR-P-3, and BR-P-16), we observe lower forecast mismatch for the experiment in Case 2. However, in some other wells (e.g., BR-P-6), the situation is the opposite, and higher forecast mismatch is observed instead. The complexity of performance comparison may be partially due to the geological structure of the reservoir model. In this benchmark case, the presence of a fault could have a substantial influence on intra-well fluid flows, which means that it could be relatively easy for the tracer to flow to wells close to the fault and relatively difficult to more distant wells. By checking the distances between the injector–producer pairs, one can obtain a clue for a possible explanation for these noticed differences in the forecast performance. For instance, BR-P-6 is located close to the south-east boundary of the reservoir, so its communication to the other parts of the reservoir may appear weaker, since tracer data from BR-I-6 is not detected. Consequently, it could lead to lower correlations between production in BR-P-6 and the petro-physical parameters in a large part of the reservoir, which makes it more difficult to improve the qualities of estimated parameters close to BR-P-6, hence leading to higher forecast mismatch.
Other than the history matching performance, we also report in Appendix B the overall wall-clock time consumed in Cases 1 and 2. As one can see there, including seismic and/or tracer data into the history matching workflow substantially increases the computational time in our research environment. Nevertheless, the overall computational time is still at a reasonable and acceptable level, meaning that our proposed workflow has the potential for real field case studies.

5. Summary and Conclusions

In this work, we investigated a joint history matching problem with multiple types of field data and proposed a coherent way to integrate production, tracer, and 4D seismic data into a history matching workflow. The workflow is demonstrated in the Brugge benchmark case using the IES-RLM algorithm, with two different experiment settings (Cases 1 and 2).
The main idea behind adding multiple types of field data is to improve the reliability of reservoir models and consequently, the forecast performance. Through two different experiments, it is shown that jointly history matching production, tracer, and 4D seismic data results in lower data mismatch for production data and lower RMSE for PERMY and PERMZ, but at the cost of slightly higher data mismatch for 4D seismic data and RMSE for PORO and PERMX. In terms of the averaged RMSE over all estimated petro-physical parameters, the reliability of the reservoir models is improved in Case 2 in comparison with Case 1.
In both cases, adopting correlation-based adaptive localization helps to maintain substantial ensemble variability even in the presence of multiple types of field data, and ensemble collapse of reservoir models is avoided. Also shown in [3], the localization scheme appears beneficial to achieve a better performance during the forecast period. By a well-to-well analysis of the data mismatch during the forecast period, it is observed that both experiments appear to have good data mismatch. Nevertheless, inter-well tracer data seem to be helpful for further reducing data mismatch, which indicates a better inter-well communication in the reservoir.
Although a better history matching performance (in terms of RMSE) is achieved with the inclusion of tracer data in Case 2, the complexity of the history matching process is increased by adding more types of field data. In the proposed workflow, we adopted a scaling procedure to adjust the relative weights among production, tracer, and 4D seismic data, which appears to work reasonably well. However, to further improve the history matching performance, it could be beneficial to develop more sophisticated methods to better balance the influence of multiple types of field data in ensemble-based history matching. Other possible lines of future research would be to consider the use of multiple tracer injectors at different reservoir locations, and to extend our proposed workflow to hydraulically or naturally fractured reservoirs [39,40]. Such investigations will be considered in our future work.

Author Contributions

Methodology, W.C.C., X.L. and K.R.P.; validation, W.C.C. and X.L.; formal analysis, W.C.C. and X.L.; investigation, W.C.C. and X.L.; data curation, W.C.C. and X.L.; writing—original draft preparation, W.C.C.; writing—review and editing, W.C.C., X.L. and K.R.P.; visualization, W.C.C.; supervision, X.L. and K.R.P. All authors have read and agreed to the published version of the manuscript.

Funding

This research was partially funded by the Research Council of Norway (project numbers: 230303 and 331644) and industry partners.

Acknowledgments

The authors would like to thank reviewers for their constructive and insightful comments and suggestions, which help to significantly improve the quality of this work. William C. Cruz acknowledges financial support from the National IOR Center of Norway (project number: 230303), which is funded by the RCN and industry partners ConocoPhillips, Aker BP, Vår Energi, Equinor, Neptune Energy, Lundin, Halliburton, Schlumberger, and Wintershall Dea. Xiaodong Luo acknowledges partial financial support from the National Centre for Sustainable Subsurface Utilization of the Norwegian Continental Shelf (NCS2030), Norway, which is funded by the Research Council of Norway (project number: 331644) and industry partners. The authors would also like to thank Schlumberger for providing academic licenses to ECLIPSE©.

Conflicts of Interest

The authors declare no conflict of interest.

Nomenclature

g Forward model
N d Number of observed data points
N e Number of reservoir models
d s i m Vector of simulated data
m Vector of reservoir model parameters
d o Vector of observed data
m true Vector of ground-truth model parameter
m ¯ Vector of mean reservoir model parameters
δ Vector of contamination noise
C d Measurement error-covariance matrix
N m Number of reservoir parameter datapoints
M Ensemble of reservoir parameters
C m Model error-covariance matrix
S m Squared-root matrix of reservoir parameters
S d Squared-root matrix of simulated data
S ˜ d Normalized squared-root matrix of simulated data
I Identity matrix
Δ d ˜ Normalized innovation
U ^ Matrix of the left singular vectors of S d
V ^ Matrix of the right singular vectors of S d
Σ ^ Matrix of the kept leading singular values
sv Kept leading singular values
N sv Number of kept singular values
K Kalman gain matrix
C Localization matrix
γ Regularization parameter
c ks Elements of C at the k-th row and the s-th column
k ks Elements of K at the k-th row and the s-th column
ρ k s Sample correlation of k-th element of the initial ensemble and s-th element of Δ d ˜
θ s Threshold value
ρ ^ s Sampling error of correlations
iIndex of iteration step
jIndex of ensemble member
K HM Effective bulk of the reservoir rock
μ HM Effective shear of the reservoir rock
ϕ c Critical porosity
μ s Gran shear modulus
v s Poisson’s ration
P eff Effective stress
C p Average number of contacts per sphere
ndegree of root
K eff Effective dry bulk modulus
μ eff Effective shear modulus
ϕ Porosity values
K s Grain bulk modulus
K f Effective fluid bulk modulus
K w Bulk modulus of water
K o Bulk modulus of oil
S w Saturation of water
S o Saturation of oil
ρ sat Saturated rock density
ρ m Mineral density
ρ w Water density
ρ o Oil density
V p P-wave velocity
V s S-wave velocity

Appendix A. Forward Simulation of 4D Seismic Data

The forward simulation of 4D seismic data goes through a few steps: reservoir simulation, petro-elastic model (PEM), and computation of seismic traces for a given time window. Starting from the initial ensemble containing the uncertain petro-physical parameters, i.e., porosity and permeability, one performs reservoir simulation to compute pore pressure and fluid saturations. Once one generates the dynamic reservoir properties, the next step is to compute the seismic elastic attributes, e.g., P-wave velocity, S-wave velocity, and formation density, by using a PEM. Finally, one has to transform the acquired seismic elastic data into another domain, e.g., amplitude-versus-angle (AVA), by using the AVA equation.
The PEM provides an important link to the reservoir model which governs the dynamics of fluid flow, and seismic elastic attributes which govern wave propagation. This connection is needed to convert pore pressure and fluid saturation into seismic elastic attributes for seismic interpretation or inversion. By far the most widely used method to establish this connection is the soft-sand model [41], in which the first step is to estimate the effective bulk ( K HM ) and shear modulus ( μ HM ) of the reservoir rock through Hertz–Minlin theory [42] as in
K HM = C p 2 ( 1 ϕ c ) 2 μ s 2 18 π 2 ( 1 v s ) 2 P eff 1 / n ,
and
μ HM = 5 4 v s 5 ( 2 v s ) 3 C p 2 ( 1 ϕ c ) 2 μ s 2 2 π 2 ( 1 v s ) 2 P eff 1 / n ,
where ϕ c , μ s , v s , and P eff represent critical porosity, gran shear modulus, Poisson’s ratio, and effective stress (e.g., lithostatic pressure minus pore pressure), respectively. In this work, the coordinate number C p , which denotes the average number of contacts per sphere, is set to 9, the degree of root n is set to 3, and ϕ c is set to 36%.
The modified Hashin–Shrikman model [43] is used to calculate the effective dry bulk modulus ( K eff ) and the effective shear modulus ( μ eff ) for porosity values ( ϕ ), as in
K eff = ϕ / ϕ c K HM + 4 3 μ HM + 1 ϕ / ϕ c K s + 4 3 μ HM 1 4 3 μ HM ,
and
μ eff = ϕ / ϕ c μ HM + μ HM 6 9 K HM + 8 μ HM K HM + 2 μ HM + 1 ϕ / ϕ c μ s + μ HM 6 9 K HM + 8 μ HM K HM + 2 μ HM 1 μ HM 6 + 9 K HM + 8 μ HM K HM + 2 μ HM ,
where K s is the grain bulk modulus.
Once one calculates the effective dry bulk modulus and the effective shear modulus for each reservoir gridblock, the next step is to compute the saturated bulk modulus and shear modulus ( K sat and μ sat , respectively) by including the saturation effect with the Gausmmann model [44], as in
K sat = K eff + 1 K eff K s 2 ϕ K f + 1 ϕ K s K eff K s 2 ,
and
μ sat = μ eff ,
where K f is the effective fluid bulk modulus, in which the two-phase fluid mixture is given by
K f = S w K w + S o K o 1 ,
where K w , K o , S w , S o , are bulk modulus of water, bulk modulus of oil, saturation of water, and saturation of oil, respectively. Further, one calculates the saturated rock density ( ρ sat ), as in [41]
ρ sat = ( 1 ϕ ) ρ m + ϕ S w ρ w + ϕ S o ρ o ,
where ρ m , ρ w , and ρ o are the mineral density, water density, and oil density, respectively. Finally, one can obtain P-wave and S-wave velocities ( V p and V s ), which can be expressed as in [41]
V p = K sat + 4 3 μ sat ρ sat ,
and
V s = μ sat ρ sat .
After the seismic elastic attributes are calculated based on the outputs of reservoir simulation, we can generate synthetic seismogram by using the Zoeppritz equation [41], in which the reflection terms between two adjacent layers are defined as a function of travel time. Then, one can obtain the desired seismic AVA data by applying a convolution between the reflectivity series and a Ricker wavelet (with a dominant frequency of 45 Hz). Here, the AVA attributes are obtained without involving any inversion process, which avoids the introduction of biases and extra uncertainties into seismic history matching later [45]. For more information regarding the procedure to generate AVA data, readers are referred to, e.g., [37].

Appendix B. Comparison of CPU Time in Different Experiments

Table A1 reports the computational cost in terms of the total wall-clock time to run the ensemble-based workflow in different experiments, with Intel(R) Core(TM) i9-10900K CPU @ 3.70GHz and 64 GB memory. For the purpose of comparison, we include the CPU time in a base case with only production data (referred to as Base case hereafter), in addition to the two cases (Cases 1 and 2) performed in this work. The results shows that including 4D seismic data (Case 1) or tracer and 4D seismic data (Case 2) considerably increase the computational time in the case studies. As one can see, in comparison to the base case, the CPU time in Cases 1 and 2 are more than doubled (+111.08 % and +120.38%, respectively). Meanwhile, in comparison to Case 1, one can also notice that the inclusion of tracer data in Case 2 slightly increases the computational cost. Overall, all case studies are finished within 30 h (108,000 s) in wall-clock time. Given that the Brugge benchmark is a field-scale reservoir model, we conclude that the computational costs in all case studies are reasonable, and that the application of the proposed workflow to other case studies at a similar scale should be affordable in general. It is important to stress that in the current work we carry out forward simulations and model updates only using a single processor. The proposed history matching workflow is, however, highly parallelizable. This means that with the aid of a high-performance-computing (HPC) facility and an optimization of our research code, one can substantially improve the computational efficiency of our proposed workflow.
Table A1. CPU time test through different experiments settings.
Table A1. CPU time test through different experiments settings.
ExperimentsBaseCase 1Case 2
CPU time (in seconds)43,30891,417 ( + 111.08 % ) 95,444 ( + 120.38 % )

References

  1. Jansen, J.D.; Brouwer, R.; Douma, S.G. Closed loop reservoir management. In Proceedings of the SPE Reservoir Simulation Symposium, Woodlands, TX, USA, 2–4 February 2009; OnePetro: Richardson, TX, USA, 2009. [Google Scholar]
  2. Lorentzen, R.J.; Bhakta, T.; Grana, D.; Luo, X.; Valestrand, R.; Nævdal, G. Simultaneous assimilation of production and seismic data: Application to the Norne field. Comput. Geosci. 2020, 24, 907–920. [Google Scholar] [CrossRef]
  3. Cruz, W.C.; Luo, X.; Petvipusit, K.R. Joint history matching of production and tracer data through an iterative ensemble smoother: A 3D field-scale case study. In Proceedings of the International Petroleum Technology Conference, Riyadh, Saudi Arabia, 21–23 February 2022; OnePetro: Richardson, TX, USA, 2022. [Google Scholar]
  4. Tayyib, D.; Al-Qasim, A.; Kokal, S.; Huseby, O. Overview of tracer applications in oil and gas industry. In Proceedings of the SPE Kuwait Oil & Gas Show and Conference, Mishref, Kuwait, 14 October 2019; Society of Petroleum Engineers: Mishref, Kuwait, 2019. [Google Scholar]
  5. Huseby, O.K.; Andersen, M.; Svorstol, I.; Dugstad, O. Improved Understanding of Reservoir Fluid Dynamics in the North Sea Snorre Field by Combining Tracers, 4D Seismic and Production Data. In Proceedings of the SPE Middle East Oil and Gas Show and Conference, Manama, Bahrain, 11 March 2007; OnePetro: Richardson, TX, USA, 2007. [Google Scholar]
  6. Allison, S.; Pope, G.A.; Sepehrnoori, K. Analysis of field tracers for reservoir description. J. Pet. Sci. Eng. 1991, 5, 173–186. [Google Scholar] [CrossRef]
  7. Datta-Gupta, A.; Lake, L.W.; Pope, G.A. Characterizing heterogeneous permeable media with spatial statistics and tracer data using sequential simulated annealing. Math. Geol. 1995, 27, 763–787. [Google Scholar] [CrossRef]
  8. Iliassov, P.A.; Datta-Gupta, A.; Vasco, D. Field-scale characterization of permeability and saturation distribution using partitioning tracer tests: The Ranger field, Texas. In Proceedings of the SPE Annual Technical Conference and Exhibition, New Orleans, Louisiana, 30 September 2001; OnePetro: Richardson, TX, USA, 2001. [Google Scholar]
  9. Valestrand, R.; Sagen, J.; Nævdal, G.; Huseby, O. The effect of including tracer data in the ensemble-Kalman-filter approach. SPE J. 2010, 15, 454–470. [Google Scholar] [CrossRef]
  10. Chen, H.; Poitzsch, M.E. Improved reservoir history matching and production optimization with tracer data. In Proceedings of the SPE Annual Technical Conference and Exhibition, Dallas, TX, USA, 24 September 2018; OnePetro: Richardson, TX, USA, 2018. [Google Scholar]
  11. Skjervheim, J.A.; Evensen, G.; Aanonsen, S.I.; Ruud, B.O.; Johansen, T.A. Incorporating 4D seismic data in reservoir simulation models using ensemble Kalman filter. In Proceedings of the SPE Annual Technical Conference and Exhibition, Dallas, TX, USA, 9–12 October 2005; OnePetro: Richardson, TX, USA, 2005. [Google Scholar]
  12. Emerick, A.A. Analysis of the performance of ensemble-based assimilation of production and seismic data. J. Pet. Sci. Eng. 2016, 139, 219–239. [Google Scholar] [CrossRef]
  13. Obidegwu, D.; Chassagne, R.; MacBeth, C. Seismic assisted history matching using binary maps. J. Nat. Gas Sci. Eng. 2017, 42, 69–84. [Google Scholar] [CrossRef]
  14. Davolio, A.; Schiozer, D.J. Probabilistic seismic history matching using binary images. J. Geophys. Eng. 2018, 15, 261–274. [Google Scholar] [CrossRef]
  15. Luo, X.; Bhakta, T.; Jakobsen, M.; Nævdal, G. Efficient big data assimilation through sparse representation: A 3D benchmark case study in petroleum engineering. PloS ONE 2018, 13, e0198586. [Google Scholar] [CrossRef]
  16. Al-Mudhafar, W.J. How is multiple-point geostatistics of lithofacies modeling assisting for fast history matching? A case study from a sand-rich fluvial depositional environment of Zubair formation in South Rumaila oil field. In Proceedings of the Offshore Technology Conference, Houston, TX, USA, 30 April–3 May 2018; OnePetro: Richardson, TX, USA, 2018. [Google Scholar]
  17. Thanh, H.V.; Sugai, Y. Integrated modelling framework for enhancement history matching in fluvial channel sandstone reservoirs. Upstream Oil Gas Technol. 2021, 6, 100027. [Google Scholar] [CrossRef]
  18. Soares, R.V.; Luo, X.; Evensen, G.; Bhakta, T. Handling Big Models and Big Data Sets in History-Matching Problems through an Adaptive Local Analysis Scheme. SPE J. 2021, 26, 973–992. [Google Scholar] [CrossRef]
  19. Evensen, G. Sequential data assimilation with a nonlinear quasi-geostrophic model using Monte Carlo methods to forecast error statistics. J. Geophys. Res. Ocean. 1994, 99, 10143–10162. [Google Scholar] [CrossRef]
  20. Nævdal, G.; Johnsen, L.M.; Aanonsen, S.I.; Vefring, E.H. Reservoir monitoring and continuous model updating using ensemble Kalman filter. In Proceedings of the SPE Annual Technical Conference and Exhibition, Denver, CO, USA, 5–8 October 2003; Society of Petroleum Engineers: Denver, CO, USA, 2003. SPE-84372-MS. [Google Scholar]
  21. Aanonsen, S.; Nævdal, G.; Oliver, D.; Reynolds, A.; Vallès, B. The Ensemble Kalman Filter in Reservoir Engineering: A Review. SPE J. 2009, 14, 393–412. [Google Scholar] [CrossRef]
  22. Skjervheim, J.A.; Evensen, G. An ensemble smoother for assisted history matching. In Proceedings of the SPE Reservoir Simulation Symposium, The Woodlands, TX, USA, 21–23 February 2011. SPE-141929-MS. [Google Scholar]
  23. Van Leeuwen, P.J.; Evensen, G. Data assimilation and inverse methods in terms of a probabilistic formulation. Mon. Wea. Rev. 1996, 124, 2898–2913. [Google Scholar] [CrossRef]
  24. Chen, Y.; Oliver, D. Levenberg-Marquardt forms of the iterative ensemble smoother for efficient history matching and uncertainty quantification. Comput. Geosci. 2013, 17, 689–703. [Google Scholar] [CrossRef]
  25. Emerick, A.A.; Reynolds, A.C. Ensemble smoother with multiple data assimilation. Comput. Geosci. 2013, 55, 3–15. [Google Scholar] [CrossRef]
  26. Luo, X.; Stordal, A.; Lorentzen, R.; Nævdal, G. Iterative ensemble smoother as an approximate solution to a regularized minimum-average-cost problem: Theory and applications. SPE J. 2015, 20, 962–982. [Google Scholar] [CrossRef]
  27. Luo, X. Novel Iterative Ensemble Smoothers Derived from A Class of Generalized Cost Functions. Comput. Geosci. 2021, 25, 1159–1189. [Google Scholar] [CrossRef]
  28. Ma, X.; Hetz, G.; Wang, X.; Bi, L.; Stern, D.; Hoda, N. A robust iterative ensemble smoother method for efficient history matching and uncertainty quantification. In Proceedings of the SPE Reservoir Simulation Conference, Montgomery, TX, USA, 20–22 February 2017; OnePetro: Richardson, TX, USA, 2017. [Google Scholar]
  29. Peters, L.; Arts, R.; Brouwer, G.; Geel, C.; Cullick, S.; Lorentzen, R.J.; Chen, Y.; Dunlop, N.; Vossepoel, F.C.; Xu, R.; et al. Results of the Brugge benchmark study for flooding optimization and history matching. SPE Reserv. Eval. Eng. 2010, 13, 391–405. [Google Scholar] [CrossRef]
  30. Engl, H.W.; Hanke, M.; Neubauer, A. Regularization of Inverse Problems; Springer: Berlin/Heidelberg, Germany, 2000. [Google Scholar]
  31. Evensen, G. Data Assimilation: The Ensemble Kalman Filter; Springer: Berlin/Heidelberg, Germany, 2009; Volume 2. [Google Scholar]
  32. Emerick, A.A.; Reynolds, A.C. History matching a field case using the ensemble Kalman filter with covariance localization. SPE Reserv. Eval. Eng. 2011, 14, 443–452. [Google Scholar] [CrossRef]
  33. Chen, Y.; Oliver, D.S. Cross-covariances and localization for EnKF in multiphase flow data assimilation. Comput. Geosci. 2010, 14, 579–601. [Google Scholar] [CrossRef]
  34. Luo, X.; Bhakta, T.; Nævdal, G. Correlation-based adaptive localization with applications to ensemble-based 4D seismic history matching. SPE J. 2018, 23, 396–427. [Google Scholar] [CrossRef]
  35. Luo, X.; Bhakta, T. Automatic and adaptive localization for ensemble-based history matching. J. Pet. Sci. Eng. 2020, 184, 106559. [Google Scholar] [CrossRef]
  36. Raniolo, S.; Dovera, L.; Cominelli, A.; Callegaro, C.; Masserano, F. History match and polymer injection optimization in a mature field using the ensemble Kalman filter. In Proceedings of the IOR 2013-17th European Symposium on Improved Oil Recovery, Saint Petersburg, Russia, 16–18 April 2013; European Association of Geoscientists & Engineers: Saint Petersburg, Russia, 2013; p. cp-342. [Google Scholar]
  37. Luo, X.; Bhakta, T.; Jakobsen, M.; Nævdal, G. An ensemble 4D-seismic history-matching framework with sparse representation based on wavelet multiresolution analysis. SPE J. 2017, 22, 985–1010. [Google Scholar] [CrossRef]
  38. Oliver, D.S.; Reynolds, A.C.; Liu, N. Inverse Theory for Petroleum Reservoir Characterization and History Matching; Cambridge University Press: Cambridge, UK, 2008. [Google Scholar]
  39. Li, J.; Pei, Y.; Jiang, H.; Zhao, L.; Li, L.; Zhou, H.; Zhao, Y.; Zhang, Z. Tracer flowback based fracture network characterization in hydraulic fracturing. In Proceedings of the Abu Dhabi International Petroleum Exhibition & Conference, Abu Dhabi, United Arab Emirates, 7–10 November 2016; OnePetro: Richardson, TX, USA, 2016. [Google Scholar]
  40. Salishchev, M.V.; Shirnen, A.A.; Khamadaliev, D.M.; Chaplygin, D.A. Logging, Well Testing and Microseismic Fracture Geometry Investigations: Mistakes, Lessons Learnt and Challenges. In Proceedings of the SPE Russian Petroleum Technology Conference, Virtual, 26 October 2020; OnePetro: Richardson, TX, USA, 2020. [Google Scholar]
  41. Mavko, G.; Mukerji, T.; Dvorkin, J. The Rock Physics Handbook: Tools for Seismic Analysis of Porous Media, 2nd ed.; Cambridge University Press: Cambridge, CA, USA, 2009. [Google Scholar] [CrossRef]
  42. Mindlin, R.D. Compliance of elastic bodies in contact. J. Appl. Mech. 1949, 16, 259–268. [Google Scholar] [CrossRef]
  43. Hashin, Z.; Shtrikman, S. A variational approach to the theory of the elastic behaviour of multiphase materials. J. Mech. Phys. Solids 1963, 11, 127–140. [Google Scholar] [CrossRef]
  44. Gassmann, F. Über die Elastizität Poröser Medien; Zürich, Inst. für Geophysik an der ETH: Zurich, Germany, 1951. [Google Scholar]
  45. Bhakta, T.; Luo, X.; Nævdal, G. Ensemble based 4D seismic history matching using a sparse representation of AVA data. In Proceedings of the 2016 SEG International Exposition and Annual Meeting, Dallas, TX, USA, 16–21 October 2016; OnePetro: Richardson, TX, USA, 2016. [Google Scholar]
Figure 1. Closed-loop reservoir development and management [1].
Figure 1. Closed-loop reservoir development and management [1].
Energies 15 06372 g001
Figure 2. Integrated history matching workflow with multiple types of field datasets.
Figure 2. Integrated history matching workflow with multiple types of field datasets.
Energies 15 06372 g002
Figure 3. A numerical model of the Brugge field with initial oil saturation and distribution of the wells.
Figure 3. A numerical model of the Brugge field with initial oil saturation and distribution of the wells.
Energies 15 06372 g003
Figure 4. Log permeability along the x-direction, with respect to realization # 2 of the initial ensemble from the Brugge field dataset.
Figure 4. Log permeability along the x-direction, with respect to realization # 2 of the initial ensemble from the Brugge field dataset.
Energies 15 06372 g004
Figure 5. Porosity with respect to realization # 2 of the initial ensemble from the Brugge field dataset.
Figure 5. Porosity with respect to realization # 2 of the initial ensemble from the Brugge field dataset.
Energies 15 06372 g005
Figure 6. Box plots of data mismatch with respect to production (left), tracer (middle), and 4D seismic (right) data at different iteration steps, where the vertical axes are in the logarithmic scale. The experiment results are represented by a colour scheme, in which the results of Case 1 correspond to red color; whereas those of Case 2 to green. Note that Case 1 does not use tracer data. Therefore, there is no corresponding result (in red) in the middle sub-figure.
Figure 6. Box plots of data mismatch with respect to production (left), tracer (middle), and 4D seismic (right) data at different iteration steps, where the vertical axes are in the logarithmic scale. The experiment results are represented by a colour scheme, in which the results of Case 1 correspond to red color; whereas those of Case 2 to green. Note that Case 1 does not use tracer data. Therefore, there is no corresponding result (in red) in the middle sub-figure.
Energies 15 06372 g006
Figure 7. Box plots of RMSE values at different iteration steps, with respect to porosity (upper left), and permeability in x-, y-, and z-directions (upper right, bottom left, and bottom right, respectively).
Figure 7. Box plots of RMSE values at different iteration steps, with respect to porosity (upper left), and permeability in x-, y-, and z-directions (upper right, bottom left, and bottom right, respectively).
Energies 15 06372 g007
Figure 8. Profiles of WOPR (first and third rows) and WWPR (second and fourth rows) from well BR-P-15 and BR-P-16, with respect to the initial ensemble (first column), and the final ensembles obtained in Case 1 (second column) and Case 2 (third column), respectively. In all sub-figures, the orange curves correspond to the production data (without measurement noise) generated by the reference model, red dots to the noisy measurements used in history matching, and blue curves to the forecast production data of respective ensembles of reservoir models.
Figure 8. Profiles of WOPR (first and third rows) and WWPR (second and fourth rows) from well BR-P-15 and BR-P-16, with respect to the initial ensemble (first column), and the final ensembles obtained in Case 1 (second column) and Case 2 (third column), respectively. In all sub-figures, the orange curves correspond to the production data (without measurement noise) generated by the reference model, red dots to the noisy measurements used in history matching, and blue curves to the forecast production data of respective ensembles of reservoir models.
Energies 15 06372 g008
Figure 9. Similar to Figure 8, but for the profiles of tracer data (WTPCW06) in wells BR-P-15 and BR-P-16, with respect to the initial ensemble (left) and the final ensembles obtained in Case 2 (right), respectively.
Figure 9. Similar to Figure 8, but for the profiles of tracer data (WTPCW06) in wells BR-P-15 and BR-P-16, with respect to the initial ensemble (left) and the final ensembles obtained in Case 2 (right), respectively.
Energies 15 06372 g009
Figure 10. Mean maps for permeability in x-, y-, and z-directions (first, second, and third rows, respectively) and porosity (last row) with respect to the initial ensemble, and the final ensembles obtained in Case 1 and Case 2, (first, second, and third columns, respectively), on Layer 2 of the reservoir model.
Figure 10. Mean maps for permeability in x-, y-, and z-directions (first, second, and third rows, respectively) and porosity (last row) with respect to the initial ensemble, and the final ensembles obtained in Case 1 and Case 2, (first, second, and third columns, respectively), on Layer 2 of the reservoir model.
Energies 15 06372 g010
Figure 11. As in Figure 10, but for STD maps instead.
Figure 11. As in Figure 10, but for STD maps instead.
Energies 15 06372 g011
Figure 12. Forecast data mismatch of WOPR and WWPR in production wells (P1–P20) with respect to the two different experiments.
Figure 12. Forecast data mismatch of WOPR and WWPR in production wells (P1–P20) with respect to the two different experiments.
Energies 15 06372 g012
Table 1. Values of data mismatch and RMSE with respect to the initial ensemble, and the final ensembles of Cases 1 and 2.
Table 1. Values of data mismatch and RMSE with respect to the initial ensemble, and the final ensembles of Cases 1 and 2.
Initial EnsembleCase 1Case 2
Production data mismatch ( 7.7861 ± 4.5592 ) × 10 3 1.1586 × 10 3 ± 1.4297 × 10 2 ( 85.1196 % ) 1.1400 × 10 3 ± 1.1465 × 10 2 ( 85.3585 % )
Tracer data mismatch ( 3.9620 ± 2.6721 ) × 10 1 6.0112 ± 9.3504 ( 84.8279 % )
4D Seismic data mismatch ( 9.2270 ± 3.4727 ) × 10 3 ( 6.3050 ± 1.6124 ) × 10 3 ( 31.6679 % ) ( 6.3960 ± 1.6680 ) × 10 3 ( 30.6817 % )
RMSE (PERMX) 1.5578 ± 0.5356 1.0284 ± 0.2160 ( 33.9838 % ) 1.0540 ± 0.2332 ( 32.3405 % )
RMSE (PERMY) 1.5524 ± 0.5315 1.0686 ± 0.2280 ( 31.1646 % ) 1.0522 ± 0.2286 ( 32.2211 % )
RMSE (PERMZ) 1.7523 ± 0.5454 1.2992 ± 0.2620 ( 25.8574 % ) 1.2412 ± 0.2276 ( 29.1674 % )
RMSE (PORO) 0.0253 ± 0.0039 0.0240 ± 0.0034 ( 5.1383 % ) 0.0242 ± 0.0035 ( 4.3478 % )
RMSE (all parameters) 1.4067 ± 0.4634 0.9866 ± 0.2011 ( 29.8642 % ) 0.9700 ± 0.1961 ( 31.0443 % )
Table 2. Parameter uncertainty reduction assessment in the different experiments.
Table 2. Parameter uncertainty reduction assessment in the different experiments.
ExperimentsSNV (PORO)SNV (PERMX)SNV (PERMY)SNV (PERMZ)
Case 1 92.43 % 69.13 % 69.79 % 75.02 %
Case 2 92.94 % 69.27 % 69.28 % 71.93 %
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Cruz, W.C.; Luo, X.; Petvipusit, K.R. Joint History Matching of Multiple Types of Field Data in a 3D Field-Scale Case Study. Energies 2022, 15, 6372. https://doi.org/10.3390/en15176372

AMA Style

Cruz WC, Luo X, Petvipusit KR. Joint History Matching of Multiple Types of Field Data in a 3D Field-Scale Case Study. Energies. 2022; 15(17):6372. https://doi.org/10.3390/en15176372

Chicago/Turabian Style

Cruz, William Chalub, Xiaodong Luo, and Kurt Rachares Petvipusit. 2022. "Joint History Matching of Multiple Types of Field Data in a 3D Field-Scale Case Study" Energies 15, no. 17: 6372. https://doi.org/10.3390/en15176372

APA Style

Cruz, W. C., Luo, X., & Petvipusit, K. R. (2022). Joint History Matching of Multiple Types of Field Data in a 3D Field-Scale Case Study. Energies, 15(17), 6372. https://doi.org/10.3390/en15176372

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop