Abstract
Developing automatic history matching (AHM) methods to replace the traditional manual history matching (MHM) approach in adjusting the permeability distribution of the reservoir simulation model has been studied by many authors. Because permeability values need to be evaluated at hundreds of thousands of grid cells in a typical reservoir simulation model, it is necessary to apply a reparameterization technique to allow the optimization algorithms to be implemented with fewer variables. In basic reparameterization techniques including zonation and pilot point methods, the calibrations are usually based solely on the production data with no systematic link to the geological and geophysical data, and therefore, the obtained permeability distribution may be not geologically consistent. Several other reparameterization techniques have attempted to preserve geological consistency by incorporating 4D seismic data; however, these techniques cannot be applied to our fractured basement reservoirs (FBRs) as they do not have 4D seismic data. Taking into account these challenges, in this study, an AHM methodology and workflow have been developed using a new reparameterization technique. This approach attempts to minimize the potential for geological nonconsistency of the calibrated results by linking the permeability to geophysical data. The proposed methodology can be applied to fields with only traditional geophysical data (3D seismic and conventional well logs). In the proposed workflow, the spatial distributions of seismic attributes and geomechanical properties were calculated and estimated from 3D seismic data and well logs, respectively. After that, a feed-forward artificial neural network (ANN) model trained by the back-propagation algorithm of the relationship between initial permeability with seismic attributes and geomechanical properties of their grid cell values is developed. Then, the calibration of the permeability distribution is performed by adjustment of the ANN model. Modification of the ANN model is performed using the simultaneous perturbation stochastic approximation (SPSA) algorithm to calibrate transmission coefficients in the ANN model to minimize the discrepancy between the simulated results and observed data. The developed methodology is applied to calibrate the permeability distribution of a simulation model of Bach Ho FBR in Vietnam. The effectiveness of the methodology is evident by comparing the historical matches with an available manually history-matched simulation model. The application shows that the proposed methodology could be considered as a suitable practical approach for adjusting the permeability distribution for FBR reservoir simulation models.
1. Introduction
Because the chance of exploring conventional reservoirs in giant fields has decreased [1, 2], unconventional petroleum reservoirs are new resources to maintain the long-term stability of oil and gas supplies. The naturally fractured reservoirs are known as important unconventional reservoirs in the world. The fractured basement reservoir (FBR) is a form of the naturally fractured reservoir. FBRs are the most common type of reservoir and contribute the most oil and gas reserves to Vietnam. The largest of FBRs in Vietnam is Bach Ho located in Cuu Long Basin [3]. Many studies have conducted overviews of geological and production characteristics for FBRs [4, 5]. These studies have shown that FBRs have low matrix porosity and permeability and that fluids are mainly stored and flow in natural fractures. Due to the complexity of the fracture system, many of the model parameters are difficult to estimate with FBRs.
One of the parameters that is the most difficult to estimate with FBRs is permeability. In general, the permeability field of a reservoir simulation model is obtained through a two-step process. In the first step (integration of static data), the permeability field in the geological model is built based on the integration of static information such as core measurement, seismic data, well logs, and stratigraphy. In the second step (integration of dynamic data), the permeability field in the simulation model (obtained by upscaling from the geological model) is calibrated during the history matching process to improve the fit between the model’s prediction and production data. In the case of FBRs, since fluids mainly flow in fractures, permeability measurements from conventional core analysis are not representative of actual effective permeability. In widely accepted practices, seismic and well log data are used to estimate permeability distribution [6, 7]. Due to the lack of representative core sample measurements, the permeability fields evaluated only from seismic and well log data are often highly uncertain and need a lot of calibration in history matching processes for simulation models of FBRs in Vietnam [8–10]. In addition, calibrations of the permeability distribution of simulation models are mainly done in a manual history matching (MHM) process and often take a long time of months with large fields like Bach Ho FBR.
Developing automatic history matching (AHM) methods to replace the traditional MHM approach in adjusting the parameters of the reservoir simulation model has been studied by many authors. Many authors have reviewed studies on AHM methods [11–13]. Permeability is a major parameter of the reservoir simulation model that is usually modified in history matching processes. However, developing suitable and reliable methods and techniques to calibrate the permeability distribution of reservoir simulation models in AHM processes faces several challenges. The first difficulty in adjusting permeability is the large number of variables because permeability values need to be evaluated at hundreds of thousands of grid cells in a typical reservoir simulation model. Therefore, it is necessary to apply a reparameterization technique to allow the optimization algorithms to be implemented with fewer variables. An up-to-date review of reparameterization techniques can be found in [14]. Examples of basic reparameterization techniques include zonation and pilot point methods. In the zonation method, the reservoir is divided into several zones, in each of which the properties are treated as uniform [15, 16]. In the pilot point method, the reservoir properties are adjusted at the selected pilot points, and the changes are extrapolated to all grid cells with kriging [17, 18]. In addition to reducing the number of variables by reparameterization techniques, another way to reduce the computational cost associated with history matching, which has been carried out in several studies, is to use proxy models (e.g., [19]). Yousefzadeh et al. presented a review of studies using proxy models in reservoir history matching [20]. Another difficulty in adjusting the permeability distribution in the reservoir simulation model is to ensure the geological consistency of the permeability field because the calibrations in the simulation model are usually based solely on the production data without a systematic link to the static data. Several reparameterization techniques have been developed to preserve geological consistency in AHM processes. For example, gradual deformation was developed to gradually change stochastic reservoir models while preserving their spatial variability [21, 22]. Other studies try to improve the geological consistency of the calibration results by incorporating 4D seismic data into the historical matching process [23–25]. Yousefzadeh and Ahmadi proposed a deep learning-based reparameterization technique that could preserve geological realism by using permeability realizations [26]. However, these AHM approaches are difficult to apply to large, complex, or non-4D seismic reservoirs.
Considering the challenges mentioned above, this study was set forth to develop a new AHM methodology and workflow to adjust the permeability field in history matching simulation models of FBRs. The new approach must be convenient for application to large, complex fields and only uses traditionally available information, including 3D seismic, well log, and production data. In this workflow, an initial feed-forward artificial neural network (ANN) model of the relationship between permeability with seismic attributes and geomechanical properties of their grid cell values is developed. The back-propagation algorithm is used to train the ANN model. Then, the calibration of the permeability distribution is performed by adjustment of the ANN model. Modification of the ANN model is performed by using the simultaneous perturbation stochastic approximation (SPSA) algorithm to calibrate transmission coefficients in the ANN model to minimize the discrepancy between the simulated results and observed data. Our AHM methodology is described in Section 2. Its application to adjust the permeability field of a simulation model of Bach Ho FBR in Vietnam is presented in Section 3.
2. Methodology
2.1. Methodology Workflow
The workflow of our AHM methodology is presented in Figure 1. The approach can start with 3D seismic, well log, and production data. The workflow can be divided into the following steps: (1)Obtain spatial distributions of initial permeability, seismic attributes, and geomechanical properties(2)Select input parameters for the ANN model from seismic attributes and geomechanical properties(3)Develop an initial ANN model of the relationship between permeability with selected seismic attributes and geomechanical properties(4)Automatic history matching the reservoir simulation model based on modifying the permeability field by calibrating the ANN model
Descriptions of the methods used to perform the above steps are given in Sections 2.2 to 2.5.
2.2. Obtain Spatial Distributions of Initial Permeability, Seismic Attributes, and Geomechanical Properties
The distribution of permeability is usually already built into the basic geological model of the reservoir. For use in a simulation model, this distribution in the geological grid is upscaled to obtain the corresponding distribution in the simulation grid. The permeability distribution on the simulation grid before history matching is used as the initial permeability distribution in the workflow.
Distributions of conventional seismic attributes can be calculated and extracted from 3D seismic data using geological modeling software. To be used in the workflow, these distributions are upscaled to the simulation grid.
To use geomechanical properties as inputs of the ANN model, we predicted distributions of vertical stress (), minimum horizontal stress (), maximum horizontal stress (), pore pressure (), Poisson’s ratio (), Young’s modulus (), uniaxial compressive strength (UCS), and internal friction coefficient (). Geomechanical properties can be calculated from well logs and core measurements along the wellbore and then interpolated to obtain their spatial distributions. In reality, laboratory testing of geomechanical properties and measurement of in situ stresses is costly and time-consuming. Furthermore, these measurement results are often available only in a few wells and at several intervals. Therefore, much research has been done to establish empirical correlations between geomechanical properties and well logs. Nowadays, it is commonly accepted that geomechanical properties may be calculated from well logs with calibration [27, 28]. Various empirical formulas and specific methods can be used to predict geomechanical parameters, depending on the data conditions and the type of reservoir. Since there are very limited measured data on the geomechanical properties of our FBRs, the prediction of geomechanical properties was performed using empirical correlations with the petrophysical well logs in the proposed methodology.
The uniaxial compressive strength (UCS) was estimated from the compressional wave transit time log (DT) using an empirical correlation obtained for granitic rock [29]: where is the P-wave velocity (), the unit of UCS is megapascals, and the unit of is meters per second.
The internal friction coefficient () was also estimated from the compressional wave sonic transit time log using the following empirical correlation [30]:
Poisson’s ratio () was calculated from the compressional and shear wave sonic transit time logs (DT and DTSM) as follows [31]: where is the S-wave velocity () and the unit of is meters per second.
Young’s modulus () was estimated from the compressional and shear wave sonic transit time logs and the density log (RHOB) using the following formula [31]:
Under conditions where density is commonly measured along with wells in the oil and gas industry, the vertical stress () at any depth is calculated by integrating the density log from the surface: where is the depth; is the depth-dependent density; the -axis is vertical and positive downward, with corresponding to the earth’s surface; and is the gravitational acceleration.
Minimum horizontal stress () was calculated from Poisson’s ratio, vertical stress, and pore pressure () according to the following equation [32]:
The maximum horizontal stress () is the most difficult to determine among the stress components. In our study, the maximum horizontal stress was calculated based on the formula of Barton et al. [33] using information from borehole image log: where is the downhole mud pressure when drilling-induced tensile fractures occur and is the wellbore breakout angle which is observed from formation microimager (FMI) logs.
An overview of methods for determining pore pressure can be found in [34]. Among them, the method of calculating pore pressure from the sonic log (which is commonly available data) was used here. This formula evolved from Eaton’s study [35] relating the pore pressure to the compressional wave sonic transit time of the following form: where is the hydrostatic pore pressure, is the measured sonic transit time by well logging, and is the normal sonic transit time obtained from the normal trend line.
2.3. Select Input Parameters for ANN Model from Seismic Attributes and Geomechanical Properties
To have an ANN model that predicts the permeability reliably and is suitable for computational resources, it is necessary to select the seismic attributes and geomechanical properties that have a high relationship with the permeability as input for the ANN model. The Pearson correlation coefficient and principal component analysis are methods commonly used for selecting input variables of ANNs [36]. Jayaweera and Aziz [37] compared the performance of these two input selection methods for an application in the development of the ANN model. Their results showed that the use of the variables selected using the principal component analysis did not contribute to improving the development of the ANN model. Meanwhile, the variables selected by the Pearson correlation coefficient were successfully used.
In this study, the selection of input parameter for the ANN model is carried out based on calculating the Pearson correlation coefficient between the seismic attributes/geomechanical properties and the initial permeability. The Pearson correlation coefficient between the dataset containing grid cell values of a potential input variable and the dataset containing corresponding grid cell permeability values is calculated using the following formula [38]: where and are the arithmetic mean of the datasets and , respectively.
2.4. Develop an Initial ANN Model of the Relationship between Permeability with Seismic Attributes and Geomechanical Properties
Feed-forward ANN is the most widely used network architecture. This type of ANN is also used in our study. A feed-forward ANN consists of one input layer, one or more hidden layers, and one output layer. Each layer contains one or more processing nodes. Every node in a layer is connected with all nodes in the previous layer. Every node in the hidden and output layers sums its weighted inputs, adds a bias, and applies an activation function to the weighted sum to get the output. In this way, the input values that pass through the network are eventually transformed into one or more output values. Two types of activation functions commonly used are the identity function and the sigmoid function, respectively, represented by the following equations:
The key step in building an ANN is a learning process in which the network is trained on a dataset consisting of input values and corresponding output values. To train a feed-forward ANN, the back-propagation algorithm is used. More detailed descriptions of the feed-forward ANN and the back-propagation algorithm can be found in many documents (e.g., [39]).
In the reparameterization technique of the proposed AHM methodology, a three-layer feed-forward ANN model representing the relationship between permeability with seismic attributes and geomechanical properties is built. Figure 2 shows the schematic diagram of a feed-forward ANN with three layers: the input layer with nodes for seismic attributes and geomechanical properties; and the output layer has one node for permeability. With this type of ANN, the final output (permeability) can be calculated using the following equation: where is the value of permeability, is the number of hidden nodes, is the number of input values, and is the th input variables (seismic attributes and geomechanical properties). The symbols , , and represent biases, weights, and activation functions, respectively, with superscript denoting the output layer, superscript denoting the hidden layer, being the index of the node in the input layer, and being the index of the nodes in the hidden layer.
Our initial ANN model is trained using the set of sample dataset that contains grid cell values of initial permeability, seismic properties, and geomechanical properties. With a selected feed-forward ANN architecture, the back-propagation algorithm is used to evaluate the weights and bias of the ANN model.
2.5. Automatic History Matching the Reservoir Simulation Model Based on Modifying Permeability Field by Calibrating the ANN Model
In AHM approaches, optimization algorithms are used to minimize an objective function that quantifies the mismatches between the simulated results and the corresponding observed values.
Typically, historical production data includes bottom-hole pressure and production rates of various fluids measured monthly at producers. Reservoir simulations are also often performed with oil rate constraints; i.e., the oil production rate is set using the corresponding historical measurements. Therefore, historical production data and simulated results may differ from each other in the values of bottom-hole pressure and water production rate of producers.
The mismatch between observation and simulation of the well bottom-hole pressure can be quantified by the root mean square error (RMSE) between the observed and simulated results of the well bottom-hole pressure of all producers and at all observation times: where is the index of the producer, is the number of the producers, is the index of the observation times, is the number of observation times of the th producer, and and are the observed and simulated bottom-hole pressure at the th observed time of the th producer.
In the same way, the mismatch between observation and simulation of the water production rate can be quantified by RMSE between the observed and simulated results of the water production rate of all the producers and at all observation times: where and are the observed and simulated water production rates of the th producers at the th observed time.
The value of the objective function in this study is calculated from the combination of RMSEs between simulation and observation of both bottom-hole pressure and water production rate as follows [40]:
The values of weighted factors and in Equation (15) are given such that the RMSE of bottom-hole pressure mismatch in the initial value of the objective function is equal to the contribution of RMSE of bottom-hole pressure.
Minimizing the objective function calculated by Equation (15) means minimizing the discrepancy between the observation and simulation of all the producers and at all observed times. The task of integrating production data to adjust the permeability distribution leads to finding the minimum of the objective function which depends on permeability distribution: where is calculated by Equation (15) and is the permeability value of grid cell .
Since the objective function represented by Equation (16) has a very large number of variables (hundreds of thousands of permeability values at the grid cells), it is necessary to use reparameterization techniques so that the optimization algorithms can be able to work efficiently. It can be seen from Equation (12) and Equation (16) that calibration of the permeability distribution can be performed by modifying the weights and bias of the selected ANN architecture. Therefore, history matching the simulation model by adjusting the permeability distribution can be performed by finding the minimum of the objective function which depends on the weights and bias of the ANN model with selected architecture:
The SPSA algorithm has been used in our study to minimize the function represented by Equation (17). Applications of SPSA algorithms to do AHM were conducted by Gao et al. [41] and Son et al. [40]. The SPSA requires only two evaluations of the objective function for the gradient approximation in each iteration of the optimization procedure, regardless of the number of parameters being optimized. A detailed description of the algorithm can be found in [42].
3. Application to Bach Ho FBR
3.1. Base Simulation Model and Initial Permeability Distribution
Bach Ho is one of the largest oil fields in Vietnam’s continental shelf located in Cuu Long Basin. The main oil reservoir of Bach Ho is a fractured basement containing fracture systems and vugs with a high degree of heterogeneity. The reservoir consists of fractured/weathered granite and granodiorite sealed by highly overpressured shales, which act both as seal and source rock. The base simulation model of Bach Ho FBR used to test the method is a model developed by the field operator using the black oil ECLIPSE 100 commercial reservoir simulation software. The size of each grid cell of the simulation model is . The simulation grid consists of grid cells in the , , and directions, respectively. Among them, 134470 grid cells are involved in flow simulation (active cell). The grid properties as the permeability of the base simulation model (initial permeability) were obtained by upscaling them from a base geological model that has also been developed by the field operator. An illustration of the initial permeability distribution is shown in Figure 3. The histogram of the initial permeability distribution of active grid cells in reservoir simulation model is presented in Figure 4.
The model has been history matched by the field operator using production data of 24 years, mainly by adjusting the permeability distribution manually. This manually history-matched model has been accepted by the field operator for forecasting future production, decision-making, and reservoir management. In our work presented here, the initial permeability distribution was modified by the proposed AHM methodology and workflow presented in Section 2. All other parameters of the model remain the same as in the manually history-matched model.
3.2. Distribution of Seismic Attributes
The distributions of nine seismic attributes (3D curvature, cosine of phase, envelope, gradient magnitude, reflection intensity, relative acoustic impedance, RMS amplitude, sweetness, and variance) were extracted from the base geological model of Bach Ho FBR and upscaling to the grid of the simulation model by the Petrel software platform. The three-dimensional distribution images of 3D curvature, envelope, and gradient magnitude attributes are shown in Figures 5–7. The three-dimensional distribution images of the other six seismic attributes are presented in Figures 8–13.
3.3. Distribution of Geomechanical Properties
To use geomechanical properties as inputs of the ANN model, distributions of vertical stress, minimum horizontal stress, maximum horizontal stress, pore pressure, Poisson’s ratio, Young’s modulus, uniaxial compressive strength, and internal friction coefficient are predicted by the method presented in Section 2.2. Firstly, computations were performed to obtain the distributions along the well of these geomechanical properties for 47 wells that have enough log data for the selected empirical correlations. Using Equations (1)–(8), wellbore distributions of the geomechanical properties were estimated for these wells. Then, the co-kriging algorithm was used to create the distribution of geomechanical properties for the entire reservoir in the geological model. The three-dimensional distribution images of vertical stress, Poisson’s ratio, and Young’s modulus are shown in Figures 14–16, respectively. The three-dimensional distribution images of the other five geomechanical properties are presented in Figures 17–21.
3.4. Selection of Input Variables for ANN Model
The ANN model was built with the output variable as permeability. The input variables for the ANN model were selected from the nine above seismic attributes and eight geomechanical properties. The choice was made on the assessment of the strength of the relationship between these potential input variables and permeability. The strengths of these relationships are evaluated based on calculating the Pearson correlation coefficient between the grid cell values of seismic attributes/geomechanical properties and the initial permeability (Equation (9)).
The calculated absolute values of the Pearson correlation coefficient between permeability and different seismic attributes and geomechanical properties are shown in Table 1. The arrangement in Table 1 is in order of the absolute value of the Pearson correlation coefficient from highest to lowest. It can be seen from Table 1 that most geomechanical properties have a higher absolute value of Pearson correlation coefficients compared to those of seismic attributes. This shows that it is perfectly reasonable to include geomechanical properties as inputs of the ANN model for the prediction of permeability.
As described in Section 2, a three-layer feed-forward back-propagation ANN is built in our workflow. The output layer of the ANN model has one node (permeability). Seismic attributes and geomechanical properties have been selected as inputs for the ANN model. The first criterion in our selection is to include attributes/properties that have a higher relationship (as quantified by the absolute value of the Pearson correlation coefficient) with permeability. However, if the correlation between two certain attributes/properties is too high, selecting both is not necessary since their effect on the output variable will be similar. From this point of view, the input parameters of the ANN model have been selected in the following order: (1)Consider including each attribute/property in the selection list based on the order of the absolute value of the Pearson correlation coefficient from high to low(2)However, if the next attribute/property has a very high absolute value of the Pearson correlation coefficient (which is quantified to be greater than 0.9) with one of the selected attributes/properties, it is not selected, but it is considered to select the next attribute/property
With the above selection criteria, studies have been performed that analyze the influence of the number of input layer nodes and the number of hidden layer nodes to determine a reasonable ANN configuration. Finally, the selected ANN model was designed with an input layer of six nodes and a hidden layer of three nodes as presented in Figure 2. Six attributes/properties were selected as inputs of the ANN model, including vertical stress, Young’s modulus, gradient magnitude, envelope, Poisson’s ratio, and 3D Curvature. These selected attributes/properties are presented in italics in Table 1.
It is noted that some attributes/properties in Table 1 were not selected even though they have a higher absolute value of the Pearson correlation coefficient than the last selected feature (3D curvature). Those attributes/properties include the following: (i)Minimum horizontal stress was not selected because its absolute value of the Pearson correlation coefficient with vertical stress (selected) is equal to 0.9890326(ii)Maximum horizontal stress was not selected because its absolute value of the Pearson correlation coefficient with vertical stress (selected) is equal to 0.9823804(iii)Uniaxial compressive strength was not selected because its absolute value of the Pearson correlation coefficient with vertical stress (selected) is equal to 0.9186830(iv)Pore pressure was not selected because their absolute value of the Pearson correlation coefficient with vertical stress (selected) is equal to 0.9243060(v)Internal friction coefficient () was not selected because its absolute value of the Pearson correlation coefficient with Young’s modulus (selected) is equal to 0.9741136(vi)RMS amplitude seismic attribute was not selected because its absolute value of the Pearson correlation coefficients with envelope seismic attribute (selected) is equal to 0.9795668(vii)Sweetness seismic attribute was not selected because its absolute value of the Pearson correlation coefficients with envelope seismic attribute (selected) is equal to 0.9966662(viii)Reflection intensity seismic attribute was not selected because its absolute value of the Pearson correlation coefficients with envelope seismic attribute (selected) is equal to 0.9707658
3.5. History Matching Result
The ANN model was built for the relationship between permeability and six selected seismic and geomechanical properties. As mentioned above, the ANN has 3 layers: the input layer consists of 6 nodes corresponding to the selected seismic attributes and geomechanical properties; the hidden layer is designed with 3 nodes; and the output layer has 1 node corresponding to the permeability. To reflect the nonlinear relationship, the sigmoid function (Equation (11)) was used to transfer information between the input layer and the hidden layer. The identity function (Equation (10)) was used to transfer information between the hidden layer and the output layer.
With the selected architecture, after training using the prepared dataset in the grid cells, the obtained ANN represents the relationship between the permeability at the grid cell () with the seismic attributes and geomechanical properties at the same grid cell () in the following form (Equation (12) when the output layer has one node):
The task of integrating production data to correct the permeability distribution leads to finding the minimum of the objective function which depends on weights and bias of the ANN model form (Equation (17) when the output layer has one node): where the value of the objective function is estimated by Equation (15).
The observed dataset was collected from the field operator. The observed dataset includes historical monthly measurements of the bottom-hole pressure and oil, gas, and water rates of 122 producers. This dataset was used in the MHM process performed by the field operator. The same historical production dataset was used in the application of our AHM methodology presented here.
As mentioned in Section 2, the SPSA algorithm was used in our computer code to minimize the objective function represented by Equation (19). Each evaluation of the objective function in the optimization algorithm requires an automatic command line run of the ECLIPSE simulator. The optimization algorithm achieves convergence after about 7 days of computation on a personal computer (PC) with 3.0 GHz CPU Intel Core i7-9700 processor. A three-dimensional image of the calibrated permeability distribution is shown in Figure 22. The calibrated permeability distribution histogram of active grid cells in reservoir simulation model is presented in Figure 23.
From Figures 3, 4, 22, and 23, it can be seen that the permeability distribution in the calibrated model is more homogeneous than the permeability distribution in the initial model. In the case of FBRs, since the fluids mainly flow in fractures, permeability measurements from conventional core analysis are not representative of the actual effective permeability. In widely accepted practices in Vietnam, a fracture “halo” modeling approach ([43, 44]) based on conventional seismic interpretation methods incorporating only the larger fault systems was employed providing permeability model [40]. This approach ignores the effect of smaller cracks and can therefore estimate the permeability distribution with higher heterogeneity than it actually is. In addition, changing the permeability values of some grid cells during manual history matching can also increase the degree of inhomogeneity of the permeability distribution. For this reason, we believe that the reduction in permeability heterogeneity after calibration is reasonable.
The simulation model using the calibrated permeability distribution (hereafter called model NEW) is compared with the manually history-matched model (hereinafter referred to as model OLD) in terms of the agreement between observed data and simulated results. The RMSEs between the observation and simulation of the bottom-hole pressure () and water production rate () of these two models are shown in Table 2. and are calculated by Equations (13) and (14), respectively. The calculations have been performed with all monthly observed times of 122 producers. Table 2 shows that RMSEs of both water production rate and bottom-hole pressure of model NEW are lower than those of model OLD. The improvement is greater in the case of water production rate in comparison with the case of bottom-hole pressure.
The values of in Table 2 indicate that model NEW is statistically better in matching with observed bottom-hole pressure compared to model OLD. The mismatch between observation and simulation of bottom-hole pressure can be quantified for each producer via RSME between the observed data and the corresponding simulated results for that producer as follows: where is RSME between observation and simulation of bottom-hole pressure of th producer, is the number of observation times of the th producer, and and are the observed and simulated bottom-hole pressure at the th observation time of the th producer.
Table 3 shows the maximum, minimum, mean, and median values of of all producers obtained by the two models. It can be seen that the minimum and median values of of model NEW are all smaller than the corresponding values of model OLD. The maximum and mean values of of the two models are equivalent. It also shows that the model NEW gives a smaller RSME of bottom-hole pressure in a larger number of producers compared to model OLD.
The values of in Table 2 indicate that model NEW is statistically better in matching with observed water production compared to model OLD. The mismatch between observation and simulation of water production rate can also be quantified for each producer via RSME between the observed data and the corresponding simulated results for that producer as follows: where is RSME between observation and simulation of water production rate ofth producer, is the number of observation times of the th producer, and and are the observed and simulated water production rate at the th observation time of the th producer.
Table 4 shows the maximum, minimum, mean, and median values of of the two models. It can be seen that the maximum, mean, and median values of of model NEW are all smaller than the corresponding values of model OLD. It also shows that model NEW gives a smaller RMSE of water production rate in a larger number of producers compared to model OLD.
When comparing the values of and obtained from the two models of all producers, it can also be seen that the number of producers having better historical matches of both water production rate and bottom-hole pressure in model NEW is superior. A total of 5 producers have RMSEs of both water production rate and bottom-hole pressure obtained from model OLD which is smaller in comparison with that of model NEW. Meanwhile, model NEW has smaller RMSEs of both water production rate and bottom-hole pressure compared with model OLD at a total of 42 producers.
Comparisons of historical matching for each producer between the two models can also be observed on time plots of water production rate and bottom-hole pressure. Due to the large number of producers, only ten plots are shown for five producers with a large cumulative liquid production in Figures 24–33. Similar figures for other producers have also been plotted. From Figures 24–33 and the comparison plots for other producers, it is also seen that the model NEW gives the agreement of both water production rate and bottom-hole pressure better in a larger number of producers in comparison with that of model OLD.
The simulated results and observed data for total field water production rate and cumulative water production are presented in Figures 34 and 35. It can be seen from these figures that model NEW gives the simulated results of the field water production that match the observed values better than the simulated results of model OLD.
Combining the results presented above, it can be said that, with both comparisons, through the RMSE values and the observations on the plots, model NEW gives simulated results that are more consistent with the observed data in comparison with model OLD. The history matching time to get model NEW is also much shorter than that of the model OLD (days versus months). This demonstrates the effectiveness of the AHM method and the reparameterization technique proposed here.
3.6. Testing of Prediction
To test the predictive ability of the models after history matching, numerical simulations were performed for the next two years and compared with the observed data updated by the field operator. Only 58 producers remained in operation during these two years. The observed dataset for testing includes historical monthly measurements of oil, gas, and water rates for all 58 producers. However, bottom-hole pressure measurements were made with only 34 producers.
Model NEW is compared in terms of the agreement between observed data and simulated results in the two-year testing period. The RMSEs between the observation and simulation of the bottom-hole pressure () and water production rate () of these two models are shown in Table 5. The calculations have been performed with all monthly observed times of all producers having measurement data in the testing period. Table 5 shows that RMSEs of both water production rate and bottom-hole pressure of model NEW are lower than those of model OLD. Similar to the history matching period, the improvement in the testing period is also greater in the case of water production rate in comparison with the case of bottom-hole pressure.
The values of in Table 5 indicate that model NEW is statistically better in testing of prediction with observed bottom-hole pressure compared to model OLD. The mismatch between observation and simulation of bottom-hole pressure production rate can also be quantified for each producer via RSME between the observed data and the corresponding simulated results for that producer (Equation (20)). Table 6 shows the maximum, minimum, mean, and median values of of producers having bottom-hole pressure in the testing period obtained by the two models. It can be seen that the minimum, mean, and median values of of model NEW are all smaller than the corresponding values of model OLD. The maximum values of of the two models are equivalent. It also shows that the model NEW gives a smaller RSME of bottom-hole pressure in a larger number of producers compared to model OLD.
The values of in Table 5 indicate that model NEW is statistically better in testing of prediction with observed water production compared to model OLD. Table 7 shows the maximum, minimum, mean, and median values of of the two models. It can be seen that the maximum, minimum, mean, and median values of of model NEW are all smaller than the corresponding values of model OLD. It also shows that model NEW gives a smaller RMSE of water production rate in a larger number of producers compared to model OLD.
Comparisons of predictive ability for each producer between the two models can also be observed on time plots of water production rate and bottom-hole pressure. Due to the large number of producers, only ten plots are shown for five producers with a large cumulative liquid production in Figures 36–45. Similar figures for other producers have also been plotted. From Figures 36–45 and the comparison plots for other producers, it is also seen that the model NEW gives the agreement of both water production rate and bottom-hole pressure better in a larger number of producers in comparison with that of model OLD.
The simulated results and observed data for total field water production rate and cumulative water production in testing period are presented in Figures 46 and 47. It can be seen from these figures that model NEW gives the simulated results of the field water production that match the observed values better than the simulated results of model OLD.
Combining the results presented above, it can be said that, with both comparisons, through the RMSEs values and the observations on the plots, model NEW gives simulated results that are more consistent with the observed data in testing period in comparison with model OLD.
4. Conclusion
This study has proposed a new automatic history matching (AHM) methodology and workflow for adjusting the permeability field of FBR’s simulation models. The proposed AHM workflow uses a reparameterization technique that allows 3D seismic and well log data to be integrated when history matching production data. In this reparameterization technique, an artificial neural network (ANN) model is developed with permeability as output. Seismic attributes and geomechanical properties that are highly correlated with permeability are selected as input for the ANN model. The grid cell values of initial permeability, selected seismic attributes, and selected geomechanical properties are used to build the ANN model. Calibration of the permeability distribution is then performed by correcting the ANN model to minimize the difference between simulated results and observed data. Correction of the ANN model is performed using an optimization algorithm to modify its weights and bias. The proposed AHM methodology has been applied to calibrate the permeability distribution for a simulation model of Bach Ho FBR of Vietnam. With the execution time tens of times shorter, the proposed AHM method gave a simulation model with better historical matching and better prediction results than the model obtained by the traditional MHM method.
The application results have shown the advantages and applicability of the proposed AHM method: (i)In the proposed method, the permeability distribution is not modified based solely on the production data but is also constrained to the static data. This approach partly preserves the geological consistency of the distribution obtained after calibration(ii)The method is very fast, easy to implement, and suitable for large reservoirs(iii)The method requires only traditional 3D seismic, well log, and production data(iv)It can also be seen that the proposed method can be used to adjust not only the permeability field but also the porosity field for FBRs. It is also expected to be suitable for other types of reservoirs under the similar data condition that core sample measurements are limited or not reliable and only traditional 3D seismic, well log, and production data are available
Abbreviations
| AHM: | Automatic history matching |
| ANN: | Artificial neural network |
| FBR: | Fractured basement reservoir |
| FWPR: | Simulated field water production rate |
| FWPRH: | Observed field water production rate in history data |
| FWPT: | Simulated field water production total |
| FWPTH: | Observed field water production total in history data |
| MHM: | Manual history matching |
| RMSE: | Root mean square error |
| SPSA: | Simultaneous perturbation stochastic approximation algorithm |
| WBHP: | Simulated well bottom-hole pressure |
| WBHPH: | Observed well bottom-hole pressure in history data |
| WWPR: | Simulated well water production rate |
| WWPRH: | Observed well water production rate in history data. |
Data Availability
The data that support the findings of this study are available from the corresponding author upon reasonable request.
Conflicts of Interest
The authors declare that there is no conflict of interest.
Acknowledgments
The research was funded by the Vietnam Academy of Science and Technology (Grant number CSCL03.01/22-23).