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
(e.g., a numerical reservoir and/or seismic simulator) which outputs a
-dimensional vector containing the simulated data
, given a
-dimensional vector of reservoir model
as the input:
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
which are obtained through the following noisy observation system:
where
stands for the ground-truth model (unknown in a real reservoir),
for the true output of the forward model, and
for a
-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
, i.e.,
.
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
whose simulated output
matches the observed data
reasonably well. In practice, when the size of the observed data is much smaller than of the reservoir model parameters (e.g.,
), 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
of reservoir models that approximately minimizes a nonlinear-least-squares (NLS) cost function, as follows:
where
i stands for the index of iteration step,
j for the index of ensemble member,
for the ensemble size, and
for the sample model error-covariance matrix with respect to the prior ensemble
, in the form of
, with the squared-root matrix
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
from its predecessor
. Here, the regularization form is known as Tikhonov regularization [
30], but one can also use other types of regularization [
27]. Meanwhile,
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 (N
d) times a factor (1 in this study).
By minimizing the cost function in Equation
3, one can obtain the following approximate model update formula:
where the square root matrix
in the model space is defined as
Similar to Equation (
5), one can define a square root matrix
in the data space as
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
and
are normalized by a square root
of
(
) to arrive at the following formulations:
In a practical implementation of the model update formula Equation (
4), numerical stability could be an issue when inverting the matrix
(especially when using 4D seismic data), since the data size
is typically much larger than the ensemble size
. 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
) by applying a Truncated Singular Value Decomposition (TVSD) to the matrix
[
24,
31]. Suppose that, through the TSVD, one obtains
where
and
are unitary matrices consisting of the kept left and right eigenvectors of
, respectively, and
is a rectangular diagonal matrix containing in its diagonal a set of retained leading singular values. The integer
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
In the current work, the regularization parameter
in Equation (
11) is chosen as
where
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
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
is significantly smaller than the sizes of reservoir model
and field data
. 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
with the Kalman gain matrix
defined as
To conduct localization in Equation (
13), one can replace the Kalman gain
by a Schur product (element-wise product) between
and a localization matrix
, as in
Here, the localization matrix
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
, 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:
where
and
stand for the elements of
and
on the
k-th row and the
s-th column, respectively,
is the
s-th element of
, and
, and
the
k-th element of the vectors
, and
, respectively.
There are different ways to compute the tapering coefficients,
, 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
dependent on the sample correlation
between the
k-th element of the initial ensemble
and the corresponding initial ensemble of the
s-th innovation element
(
). Then, they defined a threshold value
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
are independent and identically distributed (i.i.d). Due to the independence assumption, one can obtain a new ensemble
by randomly shuffling the indices
j of
. After that, one can compute the sampling errors
, which represents the correlation field between the new ensemble
and the ensemble of the
s-th innovation-data point
. Under the i.i.d assumption,
and
are independent. As a result, one can take the correlation field
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
with respect to the sampling errors
as follows [
35]:
After that, one can further compute the threshold value
by
where
is the number of elements in
. 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]
where the
k-th element
of
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)
to finally obtain the tapering coefficient
.
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 ( 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 (
) 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
where
and
denote wavelet transformation and thresholding operators, respectively, and
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
, 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
in Equation (
3) by the Schur product
, where
is a vector to be specified later, and the operator
stands for a diagonal matrix whose diagonal elements are taken from the vector
.
For production and tracer data, the elements of
are determined using the following rule:
where the operator
takes the maximum value of a given type of field data, and
stands for the type (WBHP, WOPR, WWPR or WTPCW06) of the
k-th field-data point. Note that the value
is normalized by the measurement-noise STD
associated with the data point
, 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
if
is a type of production data (WBHP, WOPR, or WWPR), and to
instead if
is a tracer data point (WTPCW06).
For seismic data, the elements of
are determined using the following rule:
where
is the wavelet coefficient belonging to the finest sub-band of
.
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:
By using Equation (
24), one can compute the scaling factor as follows:
where the overline denotes the mean value over the ensemble members,
corresponds to the mean data mismatch value of 4D seismic data, and
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
is replaced by
.
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
where
and
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
and the true one
, we use the Root Mean Squared Error (RMSE), as in
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 (
) and RMSE (
) 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.