Construction of 2D and 3D rock physics templates for quantitative prediction of physical properties of a carbonate reservoir in SW of Iran

The characterization of carbonate rocks is not straightforward, as they often experience complex diagenetic processes causing them to expose wide variations in pore types. This research aims to characterize the properties of a carbonate reservoir with a complicated porous structure through rock physics principles and tools. Two representative wells from an oil field located in SW of Iran were selected, and two-dimensional (2D) and three-dimensional (3D) rock physics templates (RPTs) were constructed by employing the appropriate rock physics models. The porosity, water saturation, and pore type are considered reservoir parameters affecting carbonate rock's elastic properties and indicating the reservoir quality. The 2D RPTs described variations in two reservoir parameters in terms of elastic properties. However, they were not able to simultaneously characterize all three reservoir parameters. The proposed 3D RPTs revealed the underlying relationship of elastic properties with pore aspect ratio, water saturation, and porosity. To validate the constructed RPTs, well logging data, scanning electron microscope images, and thin section images were utilized. The RPTs were also employed to predict the reservoir properties quantitatively, and these predictions were compared with the petrophysical data. The average errors of the predicted porosity and water saturation by 3D RPT were, respectively, 1.22% and 6.66% for well A, and 2.65% and 8.18% for well B. The 2D RPTs provided three sets of predictions for porosity and water saturation (considering three specific pore aspect ratios of 0.03, 0.1, and 0.5), all with higher average errors compared to the predictions by 3D RPT for both wells. The obtained results proved that 3D RPT could predict reservoir properties more accurately. Finally, based on the estimated values of pore aspect ratio, water saturation, and porosity using 3D RPTs, the reservoir under study was divided into distinct depth intervals, and a quality level was assigned to each interval. The introduced rock physics-based procedure for a carbonate reservoir characterization could increase the reliability in predicting the reservoir properties, enhance the ability to detect the reservoir fluid, and thereby reduce the interpretation risk.


S oil
Oil saturation S w Water saturation V p P-wave Velocity V p /V s Ratio of P-wave velocity to S-wave velocity V s S-wave Velocity α Pore aspect ratio μ calcit Shear modulus of the calcit μ clay Shear modulus of the clay μ dry Shear modulus of the dry rock μ m Shear modulus of the solid matrix μ sat Shear modulus of the fluid-saturated rock ρ calcit Density of the calcit ρ clay Density of the clay ρ f Density of the fluid of the rock ρ g Density of the gas ρ i Density of each mineral in the rock ρ m Density of the mineral matrix ρ oil Density of the oil ρ sat Density of the saturated rock ρ w Density of the water ϕ Porosity

Introduction
The rock physics template (RPT) shows the trend of changes in the reservoir properties in terms of elastic/ seismic attributes and can be utilized to characterize and evaluate the reservoir rocks. The common RPTs in rock physics are two-dimensional (2D), and the ratio of P-wave velocity to S-wave velocity (V p /V s ) is usually plotted versus P-impedance. The RPT predicts the reservoir properties, which depends on the type of the employed rock physics model to construct the RPT (Avseth and Odegaard, 2004). The RPTs can often accurately predict the reservoir properties of sandstone reservoirs. Gupta et al. (2012) analyzed a sandstone reservoir's well logs and seismic-related datasets to identify its lithology and fluid type using a specific RPT constructed based on the Kuster-Toksöz rock physics model. With the Hertz-Mindlin theory, Pan et al. (2019) constructed a three-dimensional (3D) RPT for a sandstone reservoir in North America. Using this RPT, they predicted the gas hydration saturation, porosity, and clay content. Johari and Emami Niri (2021) identified the hydrocarbon content and lithological changes within the 3D area of a sandy reservoir by establishing a 2D RPT and superimposing seismic inversion results on it. Li et al. (2022) established specific types of RPTs based on P-wave impedance and attenuation and produced RTP-based 3D cubes of seismic attributes to estimate reservoir porosity and water saturation. Although RPTs can also be used to analyze the reservoir properties of carbonate rocks, their construction for this type of lithology is a more complex task. The carbonate rocks usually have heterogeneous textures resulting from the physical, chemical, and diagenetic processes during burial history. In carbonates, the porous structure is much more complicated, and they can have different porosity types, including fracture, interparticle, intraparticle, vuggy, and cavity. The elastic attributes are dependent on the porous geometry of rocks. A complex pore network of the carbonate rocks leads to a highly scattered relation between seismic velocity and reservoir porosity . The pore geometry of the carbonate rocks is an effective parameter on elastic properties, the same as the saturation and porosity parameters. Therefore, a RPT that considers only saturation and porosity may lead to deviations in estimating the reservoir characteristics and hydrocarbon detection in the carbonate rocks (Ghasemi and Bayuk 2020). The round pores strengthen the rock structure and increase the seismic velocity. On the other hand, the penny-shaped pores weaken the rock structure and reduce the velocity. The RPTs of carbonate reservoirs should be established using the rock physics models that consider pore geometry in their equations (Pang et al. 2020). Zhao et al. (2013) considered the three parameters of fracture, pore, and cavity for a carbonate reservoir to construct a RPT with the Xu and Payne model (Xu and Payne 2009) and characterized the reservoir pore-type. Pang et al. (2019) constructed a multi-scale RPT according to the attenuation parameter in one of the carbonate reservoirs in China. They used the Biot-Rayleigh theory and predicted the porosity and water saturation. Hilman and Winardhi (2019) used the Kuster-Toksöz and Gassmann model to construct a RPT for a carbonate reservoir in the South Sumatra Basin. They suggested that the aspect ratio of the porosities and the elastic modulus of minerals should be identified for a more useful application of RPTs in carbonate reservoirs.
The 2D RPTs can only predict two of the reservoir parameters and cannot accurately interpret the properties of the carbonate reservoir. However, the 3D RPTs can predict the reservoir properties of the carbonate rocks more precisely and completely . A 3D RPT was built for a gas carbonate reservoir in China by Li and Zhang (2018). They used the Gassmann and differential effective medium (DEM) models for constructing their RPT and interpreted the reservoir's porosity, water saturation, and pore type. Pang et al. (2021) built a 3D RPT for an oil reservoir using the DEM model, which revealed the relationship between the reservoir characteristics and elastic/seismic parameters. The laboratory, well logging, and seismic data were used to calibrate the templates to estimate the various pore type volumes and the clay content. Wei et al. (2021) constructed a 3D RPT that demonstrates the effects of the temperature, differential pressure, and porosity on the V p /V s ratio, P-wave impedance, and attenuation. They established their introduced 3D RPT for a deep carbonate-type reservoir and used the Biot-Rayleigh theory implemented it to calculate the reservoir properties at seismic frequency ranges. Markus et al. (2021) generated 3D RPTs according to a rock-physics model founded on penny-shaped inclusions to predict the aspect ratio, crack, and stiff pore types of deep carbonate reservoirs of the Tarim Basin located in China.
Several types of rock physics effective medium (EM) models can be employed to characterize the elastic properties of carbonate rocks. These models account for the carbonate rocks' different pore condensations and pore geometry (Mur and Vernik 2020). Some prerequisites and assumptions exist for using the rock physics models to consider the rocks' pore geometry. Some methods (e.g., Xu and Payne 2009;Zhao et al. 2013) suppose that two parameters of saturation and porosity are known and then calculate the volumetric proportions of various pore types based on an inversion scheme. The hypothesis of some other methods (e.g., Adelinet and Le Ravalec 2015) is that carbonate rocks include only two types of pores (spherical and crack) which both are filled with water. These methods infer the spherical-type porosity and the crack density from well logs/seismic data. Ba et al. (2016) have developed a new double-porous model of clay squirt flow in which a wave-induced local fluid flow occurs between the clay aggregate micropores and the grain boundary macropores. Ba et al. (2017) presented a double double-porosity model, where the saturated pores of two different fluids overlap with pores of different compressibility. Zhang et al. (2020) set up a theoretic relation between pore structure and porosity via merging the Kuster-Toksöz and critical porosity models. This new model is used to characterize the elastic properties of the rocks with multiplex pore types. Heidari et al. (2020) investigated applying empirical and semi-empirical rock-frame models for a carbonate reservoir in one of the oilfields in southwest Iran. They used deterministic and probabilistic optimization techniques to calibrate the free parameter of each model, effective elastic moduli, and densities of the minerals and saturating fluids. Zhang et al. (2021) developed wave propagation theory in a fluid-saturated infinite porous medium containing inclusions at multiple scales. They used the DEM theory of solid composites and the Biot-Rayleigh theory for the double porosity medium.
In this research, various RPTs were constructed for the Sarvak carbonate reservoir located in southwest Iran. This reservoir has a complex pore structure and involves different types of porosities. The Voigt-Reuss-Hill average (Hill 1952), DEM analytical Zhang 2011), andGassmann (1951) models were adopted to calculate the elastic properties of the rock. As a widely used rock physics model, the Voigt-Reuss-Hill average was utilized to calculate the elastic moduli of the solid matrix. The bulk and shear moduli of the dry rock were computed using DEM analytical model. The DEM model is appropriate for high-frequency laboratory conditions. Previous studies have shown that the elastic moduli of the dry rock at the high-and low-frequency conditions are slightly different from each other, and they are roughly identical (Mavko et al. 2020). Gassmann's fluid substitution model was used to predict the changes in the bulk and shear moduli of the rock saturated with different fluids. The complexity of the pore structure in the carbonate rock may lead to the fact that they may not always be able to support the assumptions of Gassmann's relation about the connection of the pores and fluids in the porous media. For this reason, applying Gassmann's model for predicting carbonate reservoir properties is a continuous debate among researchers. Many scholars investigated the performance of employing Gassmann's model in carbonate reservoirs. They concluded that Gassmann's model could be used to determine the bulk and shear moduli of the rock saturated with fluid. Some of these researchers are: Sharma et al. (2006), Grechka (2009), Gegenhuber (2015 and Gegenhuber and Pupos (2015). Thus, Gassmann's model can be used for the elastic moduli calculations in the carbonate reservoir. After esbalishing the 2D and 3D RPTs, they are utilized to predict the reservoir properties using the elastic properties data of the rock. The available data used in this study included: well logs, thin sections, and scanning electron microscope (SEM) images. These data were further used to validate the constructed RPTs. To highlight the efficiency of the proposed RPTs, the measured and predicted reservoir properties were qualitatively and quantitatively compared with each other.

Geological setting
The oil field under study, so-called Field Y in this research, is located in southwest Iran, the westernmost part of Abadan Plain. The vast majority of Iranian oil reservoirs have been located on this plate and Dezful Embayment (Farahzadi et al. 2019;Mehrabi and Rahimpour-Bonab 2014). It is understood that Field Y is a giant long-axis anticline with a NS trend known as Arabian Trend. With two crests elongated its length, the southern elevation is more than the northern (Zadeh et al. 2019). Since Abadan Plain is placed between the southwest of Arabian Shield and northwest of Zagros fold and thrust belt (ZFTB), the geological features of this area are controlled by active and passive tectonic regimes of these two regions. Nevertheless, the dominancy of the Arabian Shield features is more evident in Abadan Plain (Khodaei et al. 2021). This field has four reservoir formations in which Sarvak limestone containing heavy oil is the most important productive reservoir of this field.
The Sarvak carbonate reservoir was the main target of this research. This formation is an important geological formation of the Bangestan group in the Zagros area. The carbonates (mainly limestones) of the Sarvak were deposited in a continental passive margin from the Late Albian to the Mid Turonian (Navidtalab et al. 2019). The Sarvak Formation is mainly composed of limestone or shaly limestone deposited in a shallow marine environment (Esrafili-Dizaji et al. 2015). The lower boundary of the Sarvak Formation is conformable with Kazhdumi Formation and is defined by overlaying clean limestones of Sarvak on Kazhdumi shales. Regionally, the upper boundary of the Sarvak Formation varies greatly. The upper boundary of Sarvak can be these items: Ilam Formation with an erosive discontinuity, Laffan Formation, and Gurpi Formation (Kalanat et al. 2021).
Various diagenetic processes affected the quality of the Sarvak Formation in the field under study and altered the primary quality of the reservoir in both positive and negative ways. These diagenetic processes include micritization, bioturbation, neomorphism, dolomitization, cementation, fracturing, dissolution, and compaction. The Sarvak reservoir specifications depend on factors such as deposition under tropical climatic conditions and early diagenesis, repeated phases of subaerial exposure due to local, regional, and global-scale tectonism, and diagenetic modification during burial (Rahimpour et al. 2012). The diagenetic and geological studies showed that different types of primary and secondary porosities were observed in the Sarvak reservoir of Field Y, such as interparticle, vuggy, intraparticle, cavity, moldic, fracture, and cracks (Malekzadeh et al. 2020).
The Sarvak formation has a good continuity in the horizontal direction, and it is homogeneous in terms of diagenetic and lithologic features . Vice versa, it is heterogeneous in the vertical direction, and the reservoir quality is better in the upper part than in the lower part. The Sarvak Formation is divided into 13 zones in terms of diversities in the diagenetic and geological facies of Field Y (Razin et al. 2010). The zones named Sarvak-1 to Sarvak-12, and Sarvak-Intra. Of course, the sequences of the Sarvak Formation may not be complete in some wells, and they may not show all of the 13 zones. The available geological studies divided Sarvak zones in the Field Y into three categories: high, medium and low reservoir quality based on the reservoir quality parameters of porosity, permeability and productivity. Sarvak-3, 5 and 8 have the high reservoir quality, Sarvak-4 and 6 have the medium reservoir quality, 2,7,9,10,11,12 and Intra have the low reservoir quality.

Available data
Among the wells in Field Y, wells A and B were selected to be studied in this research. Available data are described below.

Petrophysical well logging data
A full set of wireline logging data and petrophysical evaluations for wells A and B in the depth interval of the Sarvak formation were accessible. The well logs and petrophysical evaluations used in this research are listed in Table 1.

Thin sections
The thin sections of well A were used to study the reservoir pore types. Six hundred twenty thin sections of the whole Sarvak interval were available. These samples were analyzed with the polarization microscopy method under the plane-polarized light (PPL) and the crossed-polarized light (XPL). Investigation of thin sections under PPL and XPL lights defines the porosity and its various types. The pores are usually identified with white color in PPL light and black color in XPL light. Five images, shown in Fig. 1, are taken as representative samples, and their descriptions are reported in Table 2.

SEM studies
The scanning electron microscope (SEM) studies performed well A for 280 samples of Sarvak depth interval. These SEM images were utilized to identify the porosity types and the pore geometries. As examples of the available SEM images, four representative samples are shown in Fig. 2, and their descriptions are reported in Table 3.

Rock physics models
Three rock physics models used in this research to calculate elastic properties are the Voigt-Reuss-Hill average, DEM analytical, and Gassmann's fluid substitution models. The theory and equations of these models are briefly described below.

Voigt-Reuss-Hill average model
The Voigt-Reuss-Hill average (Hill 1952) is the arithmetic mean of Voigt's upper bound (Voigt 1890) and Reuss's lower bound (Reuss 1929). It was used to calculate the elastic moduli of the solid matrix by the following equations: where M VRH , M V, and M R are the elastic moduli computed by the Voigt-Reuss-Hill, Voigt's, and Reuss's models, respectively, M i and f i are the elastic modulus and volume percentage of each mineral component, respectively.

DEM analytical model
The dry rock's bulk and shear moduli were computed using the DEM analytical model introduced by Li and Zhang (2011). They used the equations of the DEM model (Norris et al. (1) The selected representative thin sections of well A in the PPL and XPL lights. The magnificence of thin sections A, B, C, D, and E were 30, 80, 80, 30, and 30×, respectively 1985) and provided an analytical solution to calculate the elastic moduli of the dry rock that contains the ellipsoidal pores. The elastic moduli obtained by the analytical and numerical methods were compared, and their match proved the accuracy of the analytical method. The equations of the DEM analytical model are as follows: Selected SEM images of well A. The magnitude of images A, B, C, and D was (1000, 100, 400, and 800×), respectively. The working distance of images A, B, C, and D was (21.74, 21.99, 16.91, and 20.04 mm), respectively where K dry and μ dry , respectively, are the bulk and shear moduli of the dry rock, K m and μ m are the bulk and shear moduli of the mineral matrix, respectively, and ϕ is the porosity. The constant "b" is the first derivative of P * i − Q * i phrase, constant "a" can be calculated from Eq. (12) where will be placed, and α is the pore aspect ratio, i.e., the ratio of the short axis to the long axis of the ellipsoidal pore. The required terms for the computation of P * i − Q * i and its derivatives are available in Appendix A of Li and Zhang (2011). The DEM model is favorable for the experimental conditions with high frequency.

Gassmann's model
Gassmann's model was used to calculate the bulk and shear moduli of the fluid-saturated rock (Gassmann 1951). Gassmann's model works based on these assumptions: (1) The rock is homogeneous in terms of macroscopic features.
(2) There is an interconnection between the pores. (3)  where K sat and μ sat are the bulk and shear moduli of the fluidsaturated rock, respectively, and K f is the bulk moduli of the fluid. In a fluid-saturated rock, if the fluids have a homogeneous distribution, K f is calculated using Eq. (15). If the fluids have a patchy distribution nature, K f is recommended to be calculated using Eq. (16) (Dvorkin et al. 1999): where S g , S w, and S oil are the saturations of the gas, water, and oil. K g , K w, and K oil , respectively, correspond to the gas, water, and oil bulk moduli. The following equations are related to the calculation of fluid-saturated rock density (Mavko et al. 2020): where ρ sat , ρ f, and ρ m are the density of the saturated rock, fluid of the rock, and mineral matrix, respectively, ϕ is the porosity, ρ g , ρ w , and ρ oil are, respectively, the gas, water, and oil density, f i is the volume percentage of each mineral in the rock, and ρ i is the density of each mineral in the rock. The P-wave and S-wave velocities were calculated using Eqs. where V p and V s are the P-wave and S-wave velocities, respectively.

Methodology
As illustrated in Fig. 3, a step-wise procedure is used to construct 2D and 3D RPTs. All steps of the methodology were coded in the MATLAB software version R2015a. The steps and their implementation in our case study are described hereafter: • Step 1: calculation of the elastic moduli of the solid matrix The bulk (K m ) and shear modulus (μ m ) of the solid matrix should be first computed based on the Voigt-Reuss-Hill average model ]. In our case study, calcite and clay were the two main mineral components of the Sarvak reservoir. Based on the previous studies of Heidari et al. (2020), the elastic moduli of calcite and clay of the Sarvak carbonate reservoir of field Y were chosen as K calcit = 61.41 Gpa, K clay = 18.21 Gpa, μ calcit = 22.77 Gpa, and μ clay = 5.19 Gpa. Calcite and clay volume percentages were obtained, and the bulk (K m ) and shear modulus (μ m ) of the solid matrix were then computed at each depth. For the studied depth interval, the average of the calculated matrix moduli was used for constructing the RPTs, as listed below: well A: K m = 65.13 Gpa, and μ m = 25.34 Gpa; well B: K m = 61.88 Gpa and μ m = 22.99 Gpa.
Step 2: Calculation of the dry and fluid-saturated rock elastic moduli The bulk (K dry ) and shear modulus (μ dry ) of dry rock are computed using the DEM analytical model ]. Input parameters of this model are elastic moduli of the solid matrix, porosity (ϕ), and pore aspect ratio (α). K m and μ m of the solid matrix were computed in step 1.
The bulk (K sat ) and shear modulus (μ sat ) of the saturated rock were then computed using Gassmann's fluid substitution model [Eqs. (13) and (14)]. The input parameters are K dry and μ dry (outputs of DEM analytical model), K m , porosity (ϕ), and bulk modulus of the saturating fluid (K f ). The K f was computed using Eq. (15), such that water and oil were considered two saturating fluids, and according to the equations of Batzle and Wang (1992), K w = 2.73 Gpa and K oil = 1.02 Gpa were computed for specific conditions of Sarvak reservoir in Field Y. Batzle and Wang (1992), based on empirical trends, thermodynamic equations, and several Fig. 3 The step-wise procedure to construct 2D and 3D RPTs and quantitatively predict reservoir properties using produced RPTs. These steps were codded and run in MATLAB software version R2015a published data, derived the relationships of bulk moduli, densities, velocities and viscosities of brines, oils, and gasses with temperature, pressure, and composition.

• Step 3: Calculation of the fluid-saturated rock density and seismic velocities
The fluid-saturated rock density was calculated using Eq. (17). The required input densities of the fluid (ρ f ) and solid matrix (ρ m ) were, respectively, calculated using Eqs. (18) and (19). To calculate ρ m , the volume percentages of calcite and clay were first obtained, and then their respective density values were given as ρ calcit = 2.75 g/cc and ρ clay = 2.31 g/cc (Based on the previous studies of Heidari et al. 2020). Next, for the studied depth interval of Sarvak formation, the average solid matrix density was calculated: Well A: ρ m = 2.503 g/cc, and Well B: ρ m = 2.82 g/cc. For the calculation of ρ f , ρ w = 1.03 g/cc and ρ oil = 0.77 g/cc were obtained from the equations of Batzle and Wang (1992). Finally, the P-wave and S-wave velocities were computed using Eqs. (20) and (21). The elastic properties could then be calculated using the seismic velocities and fluid-saturated rock density.

• Step 4: Construction of 2D and 3D RPTs
In this study, we constructed two types of 2D RPTs: (i) the first type is based on cross-plotting of V p versus porosity; (ii) the second type is based on the cross-plot of V p /V s versus P-wave impedance. To construct 2D RPTs, two parameters were considered variable among three reservoir properties of porosity, aspect ratio, and water saturation, and one parameter was kept constant. The changes in the elastic properties concerning the variable reservoir properties were then plotted in two dimensions.
The 3D RPT was established based on cross-plotting the P-wave impedance on the X-axis, the S-wave impedance on the Y-axis, and the density on the Z-axis. In addition, all three reservoir parameters (ϕ, α, and S w ) were set as variables, and then the changes in the elastic parameters versus the reservoir properties were plotted in three dimensions. • Step 5: Quantitative prediction of reservoir properties using RPTs Once the RPTs were constructed, the log-derived values of P-wave impedance, S-wave impedance, and density were superimposed on the constructed 2D and 3D RPTs. The 2D RPT could predict two reservoir parameters, but the 3D RPT could be employed to estimate three parameters of pore aspect ratio, water saturation, and porosity. Using the least square theory, the optimum mesh point at the closest distance to each data point was first identified. For every optimal grid point, corresponding pore aspect ratio, water saturation, and porosity values were then predicted. As the final step, the elastic properties were predicted at all depths and compared to the measured data to validate the obtained results.

Results
The step-wise procedure shown in Fig. 3 was applied during the construction of 2D and 3D RPTs for wells A and B. The qualitative and quantitative analyses of the obtained results have been described below.

2D RPTs
In this study, we constructed two types of 2D RPTs: (i) the first type is based on cross-plotting of V p versus porosity; (ii) the second type is based on the cross-plot of V p /V s versus P-wave impedance. Figure 4 demonstrates the established 2D RPTs of V p versus porosity in well A. The variable values of the reservoir properties in these RPTs are: 0.01 < α < 1 and 0 < ϕ < 20%. Figure 4a is the cross-plot of measured values of V p versus porosity for Sarvak reservoir in well A, where the projected sample points are color-coded by the water saturation percentages. The modeled curves of the pore aspect ratio have also been superimposed on the cross-plot. This RPT reveals that if the pore aspect ratio remains constant, the value of V p decreases by an increase in the porosity. If porosity remains constant, the value of V p increases by an increase in pore aspect ratio. Based on the RPT, the values of V p and porosity for each data point, its assigned pore aspect ratio could be predicted. The 0.1 < α < 0.5 were assigned for those data points with low water saturation and high productivity. To validate the predicted pore geometries, the available PPL and XPL images were utilized; five of them, as examples, were placed on the upper side of Fig. 4a. Based on the depth of thin sections, the position of each thin section on the cross-plot was determined. The 2D RPT predicted the pore geometries of thin sections based on their locations, marked in Fig. 4a. A good match is reported in Table 2 between the observed and predicted values. For example, the vuggy porosity with α = 0.2 was predicted for the depth of thin section D, which matched the observed vuggy porosity. Figure 4b is similar to Fig. 4a, but each sample point has been color-coded by the depth values. SEM images are employed to validate the obtained results; four images, as examples, are displayed in the upper part of Fig. 4b. The corresponding depth of SEM images was located on the cross plot, similar to the thin section images. Based on the position of SEM images on the cross plot, the 2D RPT predicted the pore aspect ratio of samples. The predicted pore aspect ratios matched well with the observed pore geometries, as listed in Table 3. For example, the zero porosity was predicted for the depth of SEM image "A", where an acceptable match with the observed zero porosity is visible. Figure 5 represents the 2D RPTs of well B in which the V p versus porosity with variable values of the pore aspect ratio was modeled for a reservoir fully saturated with oil. Values of 0.01 < α < 1 and 0 < ϕ < 20% were chosen as the variable reservoir properties. The measured values of V p versus porosity color-coded with gamma-ray (API) are cross-plotted in Fig. 5a. The modeled 2D RPT curves were also added to the cross-plot. The projected data points are highly scattered, indicating the pore shape's effect on the compressional velocity. For example, there are data points with velocity differences up to 950 m/s, but they show the same porosity. The values of 0.1 < α < 0.5 were predicted for data points with a low gamma-ray value and high productivity. Figure 5b is another version of this 2D RPT for well B, where sample points have been color-coded by the subzones of the Sarvak reservoir interval. This RPT revealed that well B had a complex pore geometry at different depths, and its pore aspect ratios were distributed in the range of 0.01 < α < 1. In the reservoir, cavities and vugs with high pore aspect ratios and cracks and fractures with low pore aspect ratios exist. These pore geometries affect in well A with variable pore aspect ratio curves. The well log-based data points have color codes of a water saturation (%) and b depth (m). Thin sections and SEM images are displayed with marked locations on the cross-plots the compressional velocity in this reservoir. The value of 0.1 < α < 0.5 was predicted for data points corresponding to high-quality zones, such as Sravak-3, 4, 5, 6, and 8.

2D RPTs of V p /V s versus P-wave impedance
2D RPTs of V p /V s versus P-wave impedance with variable porosity and water saturation values were also constructed for wells A and B. Three distinct rock physicsbased curves modeled for specific pore aspect ratio values (α = 0.03, α = 0.1, and α = 0.5) were also added to this cross-plot. Since the RPT could not be constructed for all values of pore aspect ratio, the values of 0.03, 0.1, and 0.5, respectively, represent fracture, interparticle and vuggy types of porosities were modeled and produced curves projected on the cross-plot. Figure 6a shows the crossplot of measured values of P-wave impedance and velocity ratio that is belonged to the Sarvak interval of well A. The porosity percentages define the color of data points in this Figure. The RPT curves for α = 0.03, α = 0.1, and α = 0.5 were identified with purple, black, and blue colors. In addition, the porosity range of 0 < ϕ < 6%, 0 < ϕ < 16%,and 0 < ϕ < 35% corresponds to curves of α = 0.03, α = 0.1, and α = 0.5, respectively. The saturation range value was chosen as 0 < S w < 100% in all curves. In the RPT of Fig. 6a, the trend for changes in elastic properties was shown in (1) The water saturation is constant, while the porosity is variable: As porosity increases, the P-wave impedance decreases, and V p /V s goes up. (2) The water saturation is variable, while porosity is constant: The P-wave impedance and V p /V s increases as water saturation increases. This RPT was constructed only for three aspect ratios, so some data points did not locate within the RPT curves. The selected thin-section images of well A are also shown in Fig. 6a. Porosity and pore aspect ratio were not predicted accurately for thin sections in Fig. 6a. Predicted values did not match the observed values in Table 2 since this RPT was not constructed for all values of pore aspect ratio. For example, the thin section C was located between the curves α = 0.03 and α = 0.1. This thin section is closer to the curve α = 0.1, so ϕ = 6.5%, α = 0.1, and S w = 100% could be predicted for that data point. The predicted values of thin section C were not consistent with the observed values in Table 2.  Figure 6b represents the constructed RPT for well B. The value ranges of the porosity were 0 < ϕ < 6.5%, 0 < ϕ < 18%, and 0 < ϕ < 40% for curves α = 0.03, α = 0.1 and α = 0.5, respectively. The saturation range value was 0 < S w < 100% in all curves. As the value of the pore aspect ratio decreases, the RPT becomes wider, and its coverage increases. The log-based data points showed that this RPT could predict porosity and water saturation relatively well, while the pore aspect ratio could not be accurately predicted.
In Fig. 7, the 2D RPT of V p /V s versus P-wave impedance was constructed for well A and B, such that the porosity and pore aspect ratio were variable, and water saturation was constant with the value of zero. Figure 7a shows the established 2D RPT for well A. In the modeled 2D RPT curves, the trends of elastic properties were plotted with Black and Blue colors. Black lines were drawn for 0.01 < α < 1. The pore aspect ratio was constant for these curves, and the porosity was variable. As porosity increased, P-wave impedance decreased, and V p /V s increased. Blue lines were drawn for ϕ = 0%, ϕ = 5%, ϕ = 15%, and ϕ = 20%, the constant parameter was porosity, and the variable parameter was pore aspect ratio. The increase in pore aspect ratio caused the increase of P-wave impedance and the decrease of V p /V s Based on the position of data points on the curves, the 2D RPT could predict the corresponding porosity and pore aspect ratio. Most data points of well A were located  A and b well B, based on the cross-plot of V p /V s ratio versus the P-wave impedance. The 2D RPT was constructed for the constant value of S w = 0% and variable porosity and pore aspect ratio values. The data points were color-coded according to their water saturation values (%). On the upper side of part a the SEM images of well A were shown within the RPT curves, and their pore aspect ratio has been distributed in the range of 0.01 < α < 1. The predicted porosity values and pore aspect ratio for displayed SEM images of well A were matched with the observed values in Table 3. For example, α = 0.23 with the vuggy geometry and ϕ = 17% were predicted for SEM image A, which was consistent with the observed value. Figure 7b is similar to Fig. 7a but is derived for well B. This RPT was constructed for a fully oil-saturated reservoir, however, most data points that did not have zero water saturation were located in the RPT. As a result, the 2D RPT could not accurately characterize the quantitative relations of elastic attributes with porosity, water saturation, and the pore aspect ratio in a complex pore system of the reservoir under study. This RPT can correctly forecast the pore aspect ratio and porosity, but fluid saturation can not. Last but not least, it should be noted that the constructed 2D RPTs could not simultaneously predict water saturation, pore aspect ratio, and porosity accurately in a carbonate reservoir with complex pore geometry.

3D RPT
As displayed in Figs. 8 and 9, the rock physics-based 3D cross-plots showing the changes in the elastic properties in the three dimensions were constructed for wells A and B. The elastic properties were shown in the constructed 3D RPTs: the P-wave impedance on the X-axis, the S-wave impedance on the Y-axis, and the density on the Z-axis. The three variable reservoir parameters were porosity, water saturation, and pore aspect ratio. Figure 8a shows a constructed 3D RPT in well A, combined with the 3D space grid of the measured values of P-wave impedance, S-wave impedance, and density. The porosity percentages define the color of data points in this Figure. The range of the reservoir properties in the 3D RPT of well A were:0 < ϕ < 22%, 0 < S w < 100%, and 0.01 < α < 0.5. By changing the reservoir properties in this RPT, their effect on the elastic properties was plotted, and thus the reservoir properties can be predicted using the elastic properties. Each data point was assigned to a 3D coordinate on the RPT, and the pore aspect ratio, water saturation, and porosity were predicted. On three-dimensional surfaces shown in Fig. 8a, blue and brown surfaces show the fully water-and oil-saturated reservoirs, respectively. The pale blue surface demonstrates the trend surface map of the Pand S-wave impedance and density with variable water saturation and porosity values and the constant value of α = 0.1. The yellow lines with the constant porosity value were plotted on the 3D surfaces for 0 < ϕ < 22%. The light blue lines with the constant pore aspect ratio were plotted on the 3D surfaces for 0.01 < α < 0.5. Some data points were located between the fully water-saturated and oil-saturated reservoir surfaces. Thus, 0 < S w < 100% were predicted for them.
The SEM image A, taken from well A, was illustrated at the top of Fig. 8a. 3D RPT predicted the properties of the SEM image based on its location, which is marked in Fig. 8a. The ϕ = 3.8%, fracture with α = 0.095 and S w = 75% were predicted for SEM image A, where the obtained results for the porosity and pore aspect ratio matched well with the observed values in Table 3. Figure 8b shows the 3D RPT of well A with − 78 degrees rotation relative to the direction shown in Fig. 8a. Some data points were located on the fully oil-saturated reservoir surface, so the S w = 0% was predicted for them. The predicted values for image B were ϕ = 11.2%, a combination of both pore and vugg with α = 0.15 and S w = 0%, and for image C were ϕ = 18.9%, vugg with α = 0.23 and S w = 0%. These predicted porosity and pore aspect ratio values were consistent with the observed values in Table 3. Figure 8c displays Fig. 8a, which was rotated for 67 degrees. Some data points were located on the fully water-saturated reservoir surface, and S w = 100% was assigned. The predicted reservoir properties by the 3D RPT can divide data into several zones with different reservoir qualities. For example, data points that their predicted reservoir properties were in the range of 0 < ϕ < 9%, 0.1 < α < 0.5, and S w = 100% should belong to the water-saturated zones in the reservoir with low porosity and poor reservoir quality. The ϕ = 0.018%, a very small crack with α = 0.05 and S w = 100%, were predicted for SEM image D displayed on the upper side of Fig. 8c. The predicted porosity and pore aspect ratio matched the observed values in Table 3 with a small error.
The 3D RPT of well B in the Sarvak reservoir is displayed in Fig. 9, where the reservoir properties range on it were:0 < ϕ < 21%, 0 < S w < 100%, and 0.01 < α < 0.5. According to the positions of the data points on the porosity and pore aspect ratio lines, their porosity and pore aspect ratio could be predicted. The predicted porosities matched well with the measured petrophysical data values, and 0.01 < α < 0.5 were predicted for the whole data points of well B. The data could be categorized based on the predicted reservoir properties. For example, the highquality zone was estimated for data points whose predicted reservoir properties were: 9 < ϕ < 21%, 0.1 < α < 0.5, and S w = 0%. The data showed that the 3D RPT could display the changes in elastic properties in terms of porosity, water saturation, and pore aspect ratio. The 3D RPT showed better accuracy in predicting the reservoir properties by considering the pore aspect ratio parameter and providing more reliable results than the 2D RPT. Finally, it was proved that the 3D RPT of a carbonate reservoir could estimate pore aspect ratio, water saturation, and porosity more accurately.

Quantitative prediction of reservoir properties using RPTs
The 2D and 3D RPTs can predict the reservoir properties quantitatively. The well log-derived data of P-wave impedance, S-wave impedance, and density of the studied wells were projected into the constructed RPTs. Using the least square theory, the optimum mesh point located at the closest distance to each data point was first identified. The corresponding values for pore aspect ratio, water saturation, and porosity for every optimal grid point could then be predicted. The whole process was coded and run in the MATLAB software version R2015a. The 3D RPT of well B Fig. 10 Comparisons of the measured and predicted reservoir properties by a 2D RPT; b 3D RPT in well A. The subzones of Sarvak are also displayed. Part a: the measured curves are marked as red, and the predicted curves according to 2D RPT for α = 0.03, α = 0.1, and α = 0.5 are purple, black, and blue, respectively. Part b: the measured and predicted curves are, respectively, drawn in red and black

Well A
The results of the quantitative predictions of the reservoir properties using constructed 2D and 3D RPTs for well A are summarized in Fig. 10. Figure 10a shows the predicted porosity and water saturation curves based on the established 2D RPT along the Sarvak depth interval of well A for three specific pore aspect ratios of α = 0.03, α = 0.1, and α = 0.5 (shown in Fig. 6a). The P-impedance and velocity ratio values for each log data point determine its position within the cross-plot of 2D RPT and corresponding values of porosity and water saturation for the associated constant pore aspect ratio (0.03 or 0.1 or 0.5) can be estimated. The average error of the predicted porosity and water saturation values based on the 2D RPT for α = 0.03 were calculated at 1.73% and 9.56%, respectively. The predicted values were in good agreement with the measured values. Only a fracture with α = 0.03 was predicted for the whole Sarvak reservoir at well A. The predicted porosity and water saturation values based on the 2D RPT for α = 0.1 also matched the observed values. The average prediction errors of the porosity and water saturation were calculated 3.48% and 8.62%, respectively. However, α = 0.1 was forecasted for all data points. The predicted porosity and water saturation mean errors were 3.14% and 7.67%, respectively, based on the 2D RPT for α = 0.5. The measured and predicted values showed a relatively good correspondence, and the vugg pore type with α = 0.5 was the only predicted pore geometry. Figure 10b displays the prediction results for the pore aspect ratio, water saturation, and porosity using the constructed 3D RPT (shown in Fig. 8) in well A. The average errors of the predicted porosity and water saturation were 1.22% and 6.66%, respectively, and the predicted values were very close to the measured values. The pore aspect ratio was also predicted for the whole depth interval. Predicting the pore aspect ratio is an important feature of 3D RPT. In addition, different types of porosities with 0.01 < α < 0.55 were predicted at each depth for each data point of well A. The interparticle pores and vuggs with 0.1 < α < 0.2 are the main types of porosity in the reservoir. Low pore aspect ratio values increase the possibility of the clay layers in the reservoir that these values belong to fracture. While the high pore aspect ratio values belonging to interparticle pores and vuggs diminish the possibility of the clay layers in the reservoir. The 3D RPT was validated with the petrophysical and geological data. The reservoir can be classified into depth intervals based on the predicted parameters according to the 3D RPT prediction, which is reported in Table 4. Figure 11 presents the quantitative predictions of reservoir properties by the 2D and 3D RPTs for well B. The 2D RPT of well B was built for α = 0.03, α = 0.1, and α = 0.5, as illustrated in Fig. 6b. The 2D RPTs for α = 0.03, α = 0.1, and α = 0.5 predicted the porosity and water saturation for data points of well B. The obtained prediction results are shown in Fig. 11a. The plotted parameters and curve colors in Fig. 11a are similar to Fig. 10a. The mean errors of the predicted porosity and water saturation according to 2D RPT for α = 0.03 are 2.73% and 10.96%, respectively. The predicted and measured values correspond approximately, but α = 0.03 was forecasted for all data points. The porosity and water saturation predictions by the 2D RPT for α = 0.1 have 2.8% and 8.60% average error values, respectively. They matched well with the measured values. However, α = 0.1 was anticipated for all data points, and the pore geometry was no longer correctly anticipated. The average errors of the porosity and water saturation predictions by the 2D RPT for α = 0.5 were 2.97% and 8.20%, respectively. A fairly good match was seen between the predicted and measured  Fig. 9, and the obtained predictions of pore aspect ratio, water saturation, and porosity from 3D RPT are shown in Fig. 11b. The heading of the axes and the colors for the curves in Fig. 11b are similar to Fig. 10b. The average errors for the predicted porosity and water saturation values were 2.65% and 8.18%, respectively. The predicted values precisely corresponded with the measured values. This estimation showed that the reservoir rock of well B has various porosities with 0.01 < α < 0.5, and the interparticle pores and vuggs with 0.1 < α < 0.25 were the main porosity in the reservoir. The conformity of the forecasted reservoir properties with the Fig. 11 Comparisons of the measured and predicted reservoir properties by a 2D RPT; b 3D RPT in well B. The subzones of Sarvak are also displayed. Part a: the measured curves are marked as red, and the predicted curves according to 2D RPT for α = 0.03, α = 0.1, and α = 0.5 are purple, black, and blue, respectively. Part b: the measured and predicted curves are, respectively, drawn in red and black Sarvak-10 0 < ϕ < 0.1 0.6 < S w < 1 0.01 < α < 0.1 Very poor 3025-3040

Well B
Sarvak-11 0.06 < ϕ < 0.18 0.4 < S w < 1 0.1 < α < 0.2 Very poor petrophysical and geological data confirmed the accuracy of the 3D RPT estimation. According to the results obtained from the 3D RPT, the reservoir could be categorized based on the depths of well B in Table 5.

Discussion
• This study has used an inclusion-based rock physics modeling approach for an oil reservoir with a limestone-type lithology with complex pore geometry. In addition to Performing rock physics analysis to diagnose reservoir rock microstructure and construct the reservoir's RPT (used for reservoir properties' estimation in undrilled areas), we have investigated the capability of constructed RPTs in reservoir properties estimation of such a heterogeneous carbonate reservoir. The following problems occurred in estimating the reservoir properties using the 2D RPTs: (1) The 2D RPT was constructed only for three aspect ratios.
(2) The 2D RPT provided three sets of predictions for the water saturation and porosity; therefore, its prediction was not unique in any depth.
(3) 2D RPT could not accurately predict the pore aspect ratio in a complex carbonate reservoir. Because the pore geometry plays an important role in the quality of the carbonate reservoir rock, the 2D RPT was not successful in this case. The 3D RPT simultaneously and accurately anticipated the pore aspect ratio, water saturation, and porosity. The predictions of 3D RPT provided precise and complete results for future studies in our reservoir with a complex pore structure. • The average errors of the predicted porosity and water saturation values based on the 2D and 3D RPTs for wells A and B are summarized in Table 6. The 2D RPTs predicted only three porosity types (each with a specific aspect ratio value) for the studied reservoir interval and provided three predictions for porosity and water saturation. Therefore, its prediction results cannot accurately interpret the properties of the studied reservoir. On the other hand, the 3D RPT quantitatively predicts the pore aspect ratio at all depths and estimates a single set of reservoir porosity and water saturation. The average errors of the predicted porosity and water saturation values are smaller; hence, the estimations based on 3D RPT can be used to interpret the reservoir properties more confidently. • The well logging data of the porosity and water saturation were available in wells A and B, and thus the predicted and measured values were compared with each other in all depths. As there was no data for the measured pore aspect ratio, the predicted pore aspect ratios could not be compared to the measured values. Consequently, they cannot be quantitatively validated. The pore aspect ratio was qualitatively validated by matching the predicted and observed values in some thin sections and SEM images.
The following approach can be suggested for the quantitative validation of the pore aspect ratio. Examine the thin sections or SEM images in all depths and measure each depth's average pore aspect ratio. These values can be plotted versus depth in the curves and compared to the predicted ones. However, this method's pore aspect ratio measurement might have high errors and is a bit time-consuming. • This study aims to use rock physics principles and tools to reduce interpretation risks/uncertainties of reservoir characteristics and help achieve a better insight into the reservoir conditions. Using the calibrated RPTs, as an interpretation/prediction tool can be very useful for describing the carbonate reservoir under study using well logs. In addition, in the case of availability of the 3D seismic inversion results (i.e., cubes of P-wave velocity, S-wave velocity, and density), they can be superimposed on the constructed RPTs, and according to the location of data points, the reservoir properties can be predicted for the exploration/development region. In a broader sense, the constructed RPTs can predict reservoir rock properties from the seismically-inverted elastic attributes in un-drilled areas of the field under study (the results of our case study will be reported in our forthcoming paper). Therefore, the planning of exploration and development well locations can be carried out with much less uncertainty.

Conclusions
• This research aimed to characterize the properties of a carbonate reservoir with a complicated porous structure by constructing the 2D and 3D RPTs. • The established 2D RPTs attribute the pore geometry's effect to the water saturation or porosity in properties prediction. Nevertheless, they could not simultaneously characterize all three reservoir parameters of pore aspect ratio (pore type), water saturation, and porosity. • The 3D RPT predicts the reservoir properties precisely and concurrently and performs better in the studied carbonate reservoir. The water saturation and porosity predictions through constructed 3D RPT matched the petrophysical data. The pore aspect ratio (pore geometry) was quantitatively predicted in all depths. Thus, the proposed 3D RPTs revealed the underlying relationship of elastic properties with pore aspect ratio, water saturation, and porosity. • Establishing a precise quantitative relation among reservoir and elastic properties leads to the following key results: (1) Increase the reliability in the prediction of the reservoir properties using well logs and seismic datas.
(2) Increase the ability to detect the reservoir rock type and fluid.
(3) Reduce the interpretation risk in exploration and development regions.