Improved Gassmann Velocity Equation to Determine the Effect of Pores Sizes in Carbonate Reservoirs Characterization

Natural oil and gas reservoirs are the main assets of petroleum exploration and production industries. Proper characterization of properties of the reservoirs and reliable estimation of their future performance is therefore of immense importance. In the future, one of the most important problems in quantitative reservoir modeling is characterization of the carbonate reservoirs. These reservoirs, as one of the major hydrocarbon settings, include heterogeneous pore spaces with unknown and irregular distributions. In this study, a carbonate reservoir in southern Iran is selected as a test bed for application of a novel characterization method. Monitoring of velocity values from sonic logs has exhibited inversion in this reservoir. We attribute this inversion to the change in pore sizes for reasons that will be explained in the paper. To obtain real values of dry rock bulk modulus as an indicator of pore sizes, assuming an identifiable model, we devised a genetic algorithm to optimize the Gassmann equation. Our results show that an appropriately designed genetic algorithm can reliably predict the values of the dry rock bulk modulus accurately. Consequently, a proposal for modification of the Gassmann equation is presented by introducing a new coefficient representing the effects of pore sizes.


Introduction
Natural oil and gas reservoirs are the main assets of petroleum exploration and production industries. Proper characterization of properties of the reservoirs and reliable estimation of their future performance is therefore of immense importance. Without suitable comprehensive methods of reservoir characterization conspicuous inaccuracies in the estimation of reservoir performance are not unlikely. These possible errors may lead to major loss of estimated revenue in hydrocarbon settings in the production process. Reservoir characterization is the process of describing various reservoir properties (mathematically and empirically) using available data in a set of reservoir models that accurately and reliably predict the reservoir performance (Jong, 2005). Typical reservoir properties have potentially a significant role in petroleum field operation and management. In a quantitative reservoir characterization study, physical quantities such as porosity, saturation, permeability, and pressure are determined numerically using available data for this purpose.
A challenging problem in quantitative reservoir modeling is characterization of the carbonate reservoirs Carbonate reservoirs, as one of the major hydrocarbon settings, include heterogeneous pore spaces with unknown and irregular distributions (from microscopic pore spaces of less than 1 mm in size to macroscopic pores of about 1 cm). Without proper determination of the distribution of pore spaces it is difficult to perform reliable characterization of the carbonate reservoirs. Many have worked on the problem of pore space detection and carbonate reservoir characterization in the past and summary of their findings is briefly presented here (see for example work by [2], [13], [28], [29], [34], Some of the developments have been an attempt in correlating the pores size with parameters such as water saturation, permeability, and porosity [1], [26], Others have studied detection of faults in a carbonate reservoir using sharp contrasts between acoustic impedances [15]. [35], developed a method for facies detection using back-propagating artificial neural networks. Other noteworthy contributions have been intelligent inversion of seismic attributes to determine carbonate facies [30], [31], using multivariable statistical procedures to determine lateral changes of porosity in a carbonate field [7], and development of relationships between porosity and seismic attributes of amplitude, phase, and frequency [17], [38], utilized amplitude variation with offset (AVO) and prestacked seismic data to obtain 162 Improved Gassmann Velocity Equation to Determine the Effect of Pores Sizes in Carbonate Reservoirs Characterization information about liquids in a carbonate formation. Most of the works done illustrate that the results cannot be completely dependable due to distribution of pore spaces and their effects on the values of reservoir properties of interest. [8], believe that pores have an effect on cementation factor in reservoirs. [16], showed that fluid saturation is an important property of hydrocarbon reservoirs that depends significantly on the pores size. [37], studied the effects of pores system on electrical conductivity of carbonates and concluded that the change in pores size can change the values of electrical conductivity and hence the values of water saturation in the reservoirs. A parameter that can reflect the characteristics of pore spaces is bulk modulus [20], defined as the value of hydrostatic pressure to the volumetric strain induced by the pressure. The implication is that bigger pores create smaller bulk modulus. Gassmann's formula can be used to convert the values of bulk modulus, porosity, shear modulus, and density to P-wave velocity [9], We propose using evolutionary algorithms in combination with Gassmann's method and material's bulk modulus to reliably assess a carbonate reservoir's characteristics. We subsequently use a carbonate reservoir in southern Iran for which the values of density, porosity, and P-wave velocity are readily available as test bed for our proposed methodology. Monitoring velocity values has exhibited inversion in this reservoir and we attribute this inversion to the change in pores sizes. The bulk modulus is optimized in a genetic algorithm, as is described in the paper, and a new coefficient that is representative of the pores sizes is introduced.

Site Geology and Petrophysical Parameters
The selected hydrocarbon reservoir contains an oil column of about 200 m thick at an approximate depth of 2850 m (from meters 2850 to 3050). Figure 1 shows a schematic view of the reservoir. The hydrocarbon level is called zone "L". Different shades of gray from white to black indicate different velocities ranging from 2500 m/s (shown in white) to 5000 m/s (in black). Drilling reports of the reservoir indicate a hydrocarbon zone of mostly pure limestones [22]. There are two wells with their log measurements available, hereinafter referred to as well No. 1 and well No. 2.The well log measurements in each well include the Density ( ρ ), the Sonic Velocity (V), the Porosity ( ϕ ), and the Water Saturation values (S W ) for each layer in the reservoir. These measurements for the oil column of well No. 1 (from meter 2850 to 3050) are shown in Figure 2.

Problem Statement and Solution Strategy
Naturally one would expect that deeper layers in a reservoir should exhibit higher velocity values due to overburden pressure of the top layers; however, when plotted the velocity values in the reservoir interval (zone "L") point to an inversion in velocity against depth from 5000 m/s at the top of the reservoir to 4000 m/s at the bottom (see Figure 3). We attribute this inversion to the pore spaces and their different sizes in the reservoir zone. Recognizing the computational power of genetic algorithms and their robustness in adaptive search and optimization we optimize the Gassmann equation using a genetic algorithm to arrive at an estimation of pores sizes and their effects.

Gassmann Theory
An interesting problem in rock physics is prediction of seismic velocities in rock layers saturated with different fluids and estimation of saturated rock velocities from dry rock, technically referred to fluid substitution problem [19], When a saturated rock layer is loaded under an increment of compressive forces, caused for instance from a passing seismic wave, an incremental pore pressure change is induced which resists the compression and results in stiffening of the rock. Gassmann introduced the low-frequency equation that relates the values of density, bulk modulus, shear modulus, and porosity to the P-wave velocity value (Gassmann, 1951): To determine velocity values, ρ and ϕ are often used from well log measurements (for example, see Figure 2). For K 0 and µ , one can refer to Castagna's five classes of limestone shown in Table 1  Geertsma suggested the following empirical equation for estimation of bulk modulus in dry rocks (Geertsma, 1961): which can be used to relate the dry bulk modulus to porosity and mineral material bulk modulus. The last parameter that needs to be determined in Gassmann equation is the effective bulk modulus of pore fluid (K fl ). Using Reuss average it is possible to obtain the equivalent bulk modulus of different phases as bellow [32].   The temperature, pressure, salinity, and composition values at the reservoir level were reported as 80℃, 35 MPa, 250000 PPM, and 18.8 API, respectively [22]. Therefore, we obtain the values of K Water and K oil as 3.8 GPa and 1.8 GPa, correspondingly.
One should note that the values of velocity in Figure 3 are from sonic logs (high frequency velocities) whereas Gassmann theory is based on low frequency data.
Fortunately, the Vertical Seismic Profiling (VSP) data is available in the area at specific intervals [24]. Accordingly, the reservoir was grouped into 6 different levels shown in Table 2. Each value in columns 4 and 5 in Table 2 is the average value in its interval obtained from Figure 2. Having the value of S W , one can determine K fl from Reuss average for each interval (Table 3). At this stage, all variables are known but K 0 , µ , and K dry . K 0 and µ can take five different values according to Table   1 (same at all six intervals), but K dry will be different in each level due to different values of ϕ (according to equation 2). The values of K 0 , µ , and K dry are determined accordingly and are shown in Table 4.
To summarize, there are five states for each level of reservoir to be evaluated in Gassman equation. Determining the velocity values at any state and comparing them to V VSP illustrates that perhaps the results from Gassmann equation can be improved. This table illustrates that in none of the states of K 0 and µ the value of velocity can be predicted precisely. This is due to inaccurately assuming K dry values as representative of pore spaces. K dry should reflect the distribution and size of the pores in a medium. Equation 2 takes into account the effects of porosity and hence involves only the effects of whole size of pores to determine K dry . Consider two cases where in one a pore has occupied the 20% of the rock compared to another case where 10 pores have each taken 2% of the rock. This is similar to levels 3 and 4 in our abovementioned example. Both of these levels have the same porosity values (see Table 2) and hence the same values of K dry ( Table 4). The only parameter that is different in these two levels is the water saturation which is higher in level 4 compared to level 3, implying higher velocities in level 4 relative to level 3. What explains lower velocities instead has to do with finer size of the pores in level 3 that cause the velocity in this level to be higher than level 4 with bigger pore spaces. To estimate the pore sizes in each level then, we use a genetic algorithm in an optimization problem to maximize Gassmann equation as a fitness function with a Gaussian mutation and a permissible range of K dry between 2 and 15 (GPa).

Genetic Algorithm
Genetic Algorithms (GAs) are optimization procedures that mimic the processes of natural evolution on machine. Algorithmically, GAs perform a stochastic search based on the mechanisms of natural selection and natural genetics [12]. They can be applied to solve a variety of optimization problems that are not well suited for standard optimization algorithms, including problems in which the objective function is discontinuous, non-differentiable, stochastic, or highly nonlinear [5], [11], [21], [27].
A genetic algorithm repeatedly modifies a population of individual solutions. At each step, it selects individuals at random from the current population to be parents and uses them to reproduce the children for the next generation. Over successive generations, the population evolves toward an optimal solution. Genetic algorithms use three major operators at each step to create the next generation of solutions from the current population: Selection rules select the individuals, called parents, which contribute to the population at the next generation. Crossover rules combine two parents to form children for the next generation. Mutation rules apply random changes to individual parents to form children.

Improved Gassmann Velocity Equation to Determine the Effect of Pores Sizes in
Carbonate Reservoirs Characterization Figure 6 illustrates a cycle of a genetic algorithm in operation. The group of chromosomes in the population will be (binary) candidates for the solution. Using an objective function the fitness values of all chromosomes are evaluated and a group of parents are selected based on their fitness values (i.e., survival of the fittest) from the current population to generate next generation of offspring. Evaluating fitness of the children with the same criterion the chromosomes in the current population are then replaced by their children. The GA cycle is repeated until a desired termination criterion is reached (convergence to a solution, for example). [33] believe that the best chromosome in the final population can become highly evolved, as a superior solution to the problem, once the process of simulated evolution is completed.
As explained before, we devise a genetic algorithm, implementing Gassmann equation as fitness function, to optimize the values of K dry as representative of pores sizes.

Implementation, Results, and Discussions
The main objective of this study is to present a formulation for better estimation of pores sizes through optimization of K dry in a given reservoir. We present our formulation using the example of Section 2. The formation in the reservoir consists of pure limestone. As such, it is natural to assume that the value of K 0 is constant over depths of interest. To determine K 0 , a genetic algorithm is used utilizing the is used as a fitness function in our genetic algorithm with the fitness limit of zero. In the GA, we used a Gaussian mutation function with the limit of generations of 200 (Gaussian function adds a random number, or mutation, chosen from a Gaussian distribution, to every entry of the parent vector [18]. Given the values of ϕ , K fl , sat ρ , and V VSP , the GA predicted the value of 63 GPa for K 0 . The best fitness function value (BFFV) after 116 generations was obtained as 0.2. Figure 7 shows the fitness function evolution over the generations.  (11) In these equations, there are two unknown variables, K dry and µ , which should be calculated. One may note that equations 6-11 are highly nonlinear and unyielding to optimality criteria and typical gradient-based optimization methods but suitable for a genetic algorithm implementation. The structure of the genetic algorithm used in this step is the same as the structure of GA used previously (to determine K 0 ). The optimal values of K dry , µ , and the corresponding fitness function for each level are obtained and shown in Table 6. Figure 8 shows the fitness variation over generations of solutions at each level.   Table 7 shows the values of estimation error for each level.  With updated values of K 0 , K dry , and µ , the cause of velocity decrease in level 4 compared to level 3 is rather obvious (recall that both levels exhibit equal porosity with a higher value of water saturation in level 4). The underlying cause, as speculated before, has to do with the reduction of K dry due to an increase in pores sizes. For other levels, Table  6 indicates: Level 1 & 2: A decrease in K dry indicates that the size of pores in level 2 is bigger than those in level 1 explaining more or less an unchanged velocity values between the two levels even though there is a decrease in porosity and an increase in the water saturation.
Level 2 & 3: A large increase in K dry imply smaller pores sizes in level 3 compared to level 2 explaining why the velocity has increased between these two layers even though there is an increase in porosity and a decrease in water saturation levels.
Level 4 & 5: An increase in K dry imply smaller pores sizes in level 5 compared to level 4 explaining why the velocity has decreased slightly between these two layers compared to the 5% increase in the porosity and larger reduction of water saturation levels.
Level 5 & 6:A large reduction of K dry may indicate that the size of pores in level 6 is bigger than those in level 5 explaining reduction in velocity values between the two levels even though there is a decrease in porosity and 10% increase in the water saturation levels.
To further validate our results from the genetic algorithm solution we used Dunham's classification of the carbonate rocks (Table 8). Dunham classification includes six classes of rocks according to their sort and crystallization. These classes are Mudstone, Wackestone, Packstone, Grainstone, Boundstone and Crystalline [6]. Mudstone, on the top of the list, is considered a rock type with the smallest pores of microscopic scale. Crystalline, at the bottom of the list, is a rock type with the largest or macroscopic scale pores [25], [36]. Table 8 reinforces our genetic algorithm's results. One may ask if the values of µ can also be representative of the pore spaces. In Table 6 the value of µ has been reduced between levels 4 and 5. To justify the velocity changes between these two zones the value of µ has to be increased.
Therefore, in our study, the only parameter which can be assumed as a representative of pores sizes is K dry .. The K dry value in Gassmann equation (from Geertsma) reflecting the effect of porosity, can be replaced by the new value of K dry from GA. The updated value of K dry from GA contains the effects of porosity and pores sizes together. Comparing these two values of K dry in Table 9, it is interesting to note that the values obtained from GA are less than those from Geertsma's equation. This is attributed to the distribution and the size of the pores. Therefore, we propose a modified Gassmann equation as:  (13) where, K Gdry is the K dry value from Geertsma equation and α is a coefficient of pores sizes. The value of α is between 0 and 1 (0 < α < 1). A near zero value for α indicates that porosity of the rock is near macroscopic while close to one values indicate microscopic porosities in the rock.
Another interesting observation in our results is about the values of water saturation in levels with small values of α (levels of macroscopic large pore spaces). In our test bed, levels 2, 4, and 6 had a higher value of water saturation compared to the other 3 levels (as seen in Table 2). The hypothesis which is propounded here is the capability of saline water in dissolving carbonate rocks. Since hydrocarbon reservoirs are formed during a very long time water can act as a dissolver during their formation. Therefore, where saturation of saline water is higher, the size of pores can be expected to be larger. This hypothesis opens a new horizon in making better predictions of reservoir characteristics especially when hydrocarbon saturation is in play.

Application of the Proposed Methodology to a Different Case
To test the applicability of our proposed methodology to another realistic case, we gathered required data from well log measurements in well No. 2. The oil column was located at the depth of 2744 m with a thickness of about 260 m (from meters 2744 to 3010). Drilling reports indicated hydrocarbon zones of mostly pure limestone similar to well No. 1 [22]. In four levels of the oil column considered in well No. 2 the changes in porosity and water saturation do not show a clear relationship with the corresponding velocity changes (Table  10). In Table 10, there is no considerable change in porosity between levels 1 and 2. A slight increase in the water saturation levels warrants an increase in velocity at level 2. However, the velocity is decreased in this level. Again, bigger sizes of the pores in this level compared to level 1 can be the underlying cause. Between levels 3 and 4, there is some decrease in the porosity with large increases in the water saturation levels. The sum effect should point to increased velocity in level 4. Although the velocity in level 4 is higher than in level 3, it is not as much as one would have expected. We believe that the pores with large sizes can be responsible for this discrepancy too. Dunham's classification of the well confirms our prediction again. Table 11 shows the results of Dunham classification.

Concluding Remarks
Unknown heterogeneous distribution of pore spaces in carbonate hydrocarbon reservoirs exerts difficulties in performing reservoir characterization. Reservoir velocity, as an important indicator of the reservoir properties, is susceptible to problems caused by the pore spaces distribution. Consequently, methods of optimization such as the one proposed in this paper, may be deemed necessary in estimation of the reservoir characteristics. We successfully implemented and tested a machine learning technique (a genetic algorithm) that considers unknown heterogeneous distribution of the pore system in an optimization problem (foreseeing of the size of pores in the hydrocarbon column in a well). Our approach uses the Gassmann equation as a fitness function for the genetic algorithm. The optimization problem seeks the zero value of the objective function and strives to obtain a desirable output which is, in our case, the real values of the dry rock bulk modulus (K dry ) as an indicator of the pores size.
Two wells were considered to test our methodology. Both wells exhibited the power of genetic algorithms in optimizing the Gassmann equation reliably. Correlation between the pores sizes obtained from the genetic algorithm and those obtained from Dunham classification were observed. Using the results of this work, Gassmann equation was modified by introducing a new coefficient. The original Gassmann equation only considers the effect of pores and not their sizes. In our proposed form of the Gassmann equation, both of these effects are considered.
An interesting aspect of this study is the observed relationships between the size of the pores and the water saturation values. Results of both wells indicated that in levels with bigger pores sizes the values of water saturation were also higher. The hypothesis presented was the capability of saline water in dissolving carbonate rocks. This hypothesis can open a new horizon in making better prediction of the reservoirs characteristics especially when hydrocarbon saturation is in play. It is anticipated that a future study which is underway may focus on: 1) Performing the modified Gassmann equation for different sizes of pores for constant porosity values 170 Improved Gassmann Velocity Equation to Determine the Effect of Pores Sizes in Carbonate Reservoirs Characterization and evaluating the effects of pores on the velocity values as an important seismic attribute. Consequently, performing seismic forward modeling and extracting other seismic attributes relating to the pores sizes. 2) Given a set of seismic attributes specific to a given pores size setting, use the results of this study to produce fuzzy rules that relate the attributes of pores, porosity, and velocity to the values of water saturation and hydrocarbon saturation.