Elastic Properties of Pannonian Basin Limestone under Different Saturation Conditions

: Understanding elastic properties of reservoir rocks is essential for seismic modeling under different saturation conditions as well as lithology discrimination. Experiments on elastic properties of limestones are signiﬁcantly less published compared to siliciclastic sedimentary rocks. The current study presents the results of laboratory measurements on Pannonian Basin limestone cores. The research was carried out for the ﬁrst time for a hydrocarbon reservoir in the Bjelovar Depression, located in the southern part of the Pannonian Basin. Ultrasonic velocity measurements and determination of dynamic elastic properties were performed on limestone plugs, in dry and saturated condition under different conﬁning pressure steps. Based on the results obtained in laboratory conditions, an empirical relationship between shear wave velocity ( Vs ) and compressional wave velocity ( Vp ) has been deﬁned. The saturated samples show an effect of shear modulus weakening, while three samples have a shear modulus strengthening effect. Two models were used in the interpretation of the measured data, the Kuster and Toksöz and the Xu-Payne model. The results show that the Xu-Payne model describes the data well and the dominant pore type system in the limestone samples can been identiﬁed. The interpretation revealed an interparticle pore system with a fraction of microcracks from 20% to 35%. The results have helped to understand the elastic properties of limestones from the southern part of the Pannonian Basin, which are necessary for any process of reservoir characterization, such as porosity distribution and permeability variation. shear modulus versus the difference between measured and Gassmann-predicted bulk modulus.


Introduction
The elastic property of rock is useful information for oil or gas reservoir characterization, especially for carbonate reservoirs. Elastic parameters are essential for the reservoir seismic response modeling under different saturation conditions, also for lithology discrimination. Carbonates contain 60% of the world's oil reserves [1], but there are only a few experiments dealing with elastic properties of carbonates compared to siliciclastic sedimentary rocks. This paper presents the results of laboratory measurements of ultrasonic wave velocities on limestone cores, which were performed for the first time for a hydrocarbon reservoir located in the Pannonian Basin, Croatia. The ultrasonic velocity measurements were performed on the Bjelovar Depression gas reservoir core samples. Laboratory measurements of porosity and velocity under different confining pressures were analyzed.
According to previously published papers, e.g., [2][3][4][5] laboratory measurements have shown that velocity in water-saturated carbonates is influenced by porosity and pore shape. These studies have shown that velocity generally decreases with increasing porosity. Moreover, large velocity variability in carbonates has been observed. Rafavich et al., (1984) [2] analyzed the ultrasonic velocities of P-and S-waves in carbonate samples from four wells with different porosity and mineralogy, such as limestones, dolomites and anhydrites. The tests were performed for 98 dry and water-saturated samples under effective stress up to 40 MPa. They concluded that porosity and density are the main factors affecting velocity. According to Anselmetti and Eberli, (1993) [6], who studied the P-and S-wave velocities in carbonate samples from different reservoirs, the influence of mineral composition in carbonates is minimal. Moreover, rocks with primary porosity have lower velocity values compared to rocks with secondary porosity.
Bakhorji, (2010) [7] compared the velocity results in carbonates with the theoretical Gassmann and Biot models [8,9] and concluded that the Gassmann model consistently overestimates the velocity values of saturated samples at low effective pressure and the samples at high effective velocities are approximately equal. The Biot model overestimates the values of saturated samples in most of the analyzed samples.
The saturation state of the samples affects the elastic parameters, and accordingly the bulk modulus values of the saturated samples are higher than those of the dry samples. The shear modulus of water is zero and theoretically the shear modulus of the saturated sample should remain unchanged. However, many published papers, e.g., [5,6,10] report that the shear modulus values of saturated samples are lower than the shear moduli of dry samples, and this phenomenon is called the "shear weakening effect". According to Baechle et al., (2009) [5] and Khazanehdari and Sothcott, (2003) [11] the causes of a shear weakening effect are viscosity, reduction of free surface energy on mineral grains, and velocity dispersion due to local flow.
The present study follows and complements the research performed so far on cores from different basins around the world and presents the first measurements of ultrasonic wave velocities on limestone cores under different effective stress values in the southern part of the Pannonian Basin. Based on the results of laboratory measurements, an empirical relationship between shear wave velocities (Vs) and compressional wave velocities (Vp) has been defined. Moreover, an analysis of the elastic parameters of the limestone cores is presented.

Geological Settings of the Bjelovar Depression
The research area belongs to the Bjelovar Depression, located in the southern part of the Pannonian Basin ( Figure 1). The geological and structural development of the research area is related to the geological development of the Pannonian Basin.
The opening of the Pannonian Basin starts in the Lower Miocene during the continental collision and subduction of the European and African plate. Extension of the Pannonian Basin is enabled by two geodynamic processes: subduction of the European plate under the Carpathians, followed by slab retreat and collision of the Adriatic microplate and the European plate [12][13][14]. The syn-rift phase presented with the formation of half-graben lasted from the Ottnangian to the Middle Badenian [15].
Extension tectonics and the opening of a series of listric faults formed four elongated basins which represented the main depositional area, such as the Drava Depression and the Bjelovar Depression, which together with the Požega Valley formed a single half-graben, as well as the Sava Depression and the Karlovac Depression [16]. The extensional tectonics continued until the Middle Badenian, when regional marine transgression occurred [17,18].
During the Pliocene, new and older fault systems were activated and reactivated in the Bjelovar Depression area, resulting in the uplift of the Bilogora Mt and partly the Slavonian Mt. (Figure 1). Due to the convergence of the Adriatic microplate during the Quaternary, with counterclockwise rotation [19][20][21], tectonic activity continues.
Unlike the other parts of the Drava Depression and the Sava Depression, the Bjelovar Depression was not on the main path of material distribution during the post-extension phase, so the sediment supply was significantly lower. The entire depositional area of the Bjelovar Depression was more stagnant than the Drava Depression, and accordingly the thickness of Neogene-Quaternary deposits is significantly smaller than that of the Drava Depression [22]. The research area consists of Pliocene and Quaternary clastic sediments which transgressively overlay on Neogene limestone basement. The limestone was deposited on a carbonate platform, but it is not a classical platform carbonate. The depositional area included shallow water with higher water energy and higher biogenic material input as well as the area with lower water energy and constant influence of the open sea. Barrier reefs in the shallower part do not represent a significant barrier and formation of a lagoon, but present shelter from tidal currents, rough swell of waves and occasional storms. Around such reefs, water circulation is much easier, so facies features are often mixed. Limestone is defined as recrystallized grainstone with included bioclastic detritus presented with shallow marine fossils. A stylolite and irregularly oriented cracks are evident on the samples. According to SEM (Scanning Electron Microscopy) carbonate component consist exclusively of limestone in the form of rhombohedral calcite crystals and the stylolites are mainly filled with clay minerals such as sericite, illite and chlorite [23]. thickness of Neogene-Quaternary deposits is significantly smaller than that of the Drava Depression [22]. The research area consists of Pliocene and Quaternary clastic sediments which transgressively overlay on Neogene limestone basement. The limestone was deposited on a carbonate platform, but it is not a classical platform carbonate. The depositional area included shallow water with higher water energy and higher biogenic material input as well as the area with lower water energy and constant influence of the open sea. Barrier reefs in the shallower part do not represent a significant barrier and formation of a lagoon, but present shelter from tidal currents, rough swell of waves and occasional storms. Around such reefs, water circulation is much easier, so facies features are often mixed. Limestone is defined as recrystallized grainstone with included bioclastic detritus presented with shallow marine fossils. A stylolite and irregularly oriented cracks are evident on the samples. According to SEM (Scanning Electron Microscopy) carbonate component consist exclusively of limestone in the form of rhombohedral calcite crystals and the stylolites are mainly filled with clay minerals such as sericite, illite and chlorite [23].

Methodology
Ultrasonic velocity measurement is based on measuring the travel time of an ultrasonic wave passing through a sample. A typical ultrasonic measuring system consists of a wave source, a transducer, a receiver, and a recording and display device ( Figure 2). The wave source generates a high-voltage electrical pulse that is converted into compressional wave or polarized shear waves by the piezoelectric crystals in the transducer. The piezoelectric crystals in the receivers are used to convert acoustic energy into voltage. Ultrasonic velocity measurements were performed with a device that enables velocity measurements and determination of dynamic elastic parameters at triaxial loading of samples up to 70 MPa, using two piezoelectric crystals with a frequency of 1 MHz that generate Pwaves and polarized S-waves.
High-frequency ultrasonic waves are generated at one end of the sample. The wave propagating through the sample is detected by a receiver on the other side of the sample. The wave is converted into an electrical signal at the opposite end of the sample by a transformer and recorded with a digital oscilloscope.
Prior to velocity measurements in the laboratory, it is necessary to define the correction parameters that should be used to accurately calculate the dynamic elastic parameters

Methodology
Ultrasonic velocity measurement is based on measuring the travel time of an ultrasonic wave passing through a sample. A typical ultrasonic measuring system consists of a wave source, a transducer, a receiver, and a recording and display device ( Figure 2). The wave source generates a high-voltage electrical pulse that is converted into compressional wave or polarized shear waves by the piezoelectric crystals in the transducer. The piezoelectric crystals in the receivers are used to convert acoustic energy into voltage. Ultrasonic velocity measurements were performed with a device that enables velocity measurements and determination of dynamic elastic parameters at triaxial loading of samples up to 70 MPa, using two piezoelectric crystals with a frequency of 1 MHz that generate P-waves and polarized S-waves.
High-frequency ultrasonic waves are generated at one end of the sample. The wave propagating through the sample is detected by a receiver on the other side of the sample. The wave is converted into an electrical signal at the opposite end of the sample by a transformer and recorded with a digital oscilloscope.
Prior to velocity measurements in the laboratory, it is necessary to define the correction parameters that should be used to accurately calculate the dynamic elastic parameters of The research was performed on 8 samples of core plugs with a diameter of 1" (2.54 cm) from two boreholes which penetrated the Bjelovar Depression limestone gas reservoir ( Figure 3). The depth range of the investigated core samples is from 805.24 to 974.11 m below the surface which is taken into account during the definition of the applied values of effective pressure.
Prior to the measurement the samples were subjected to a cleaning process which involved organic solvent extraction and drying. Following the sample preparation [25], the effective porosity of the plugs was determined with the Helium expansion method at different confining stress values, namely 5, 7, 9, 11, and 13 MPa.
The ultrasonic velocity measurements were performed on the dry samples in the first place, then the plugs were vacuum saturated with 10 g/L NaCl solution which corresponds to the reservoir brine salinity. The tests were performed at multiple isostatic confining pressure values in both cases, and the applied stresses were identical with the ones used in the porosity determination. A pore pressure was not applied during measurements of the saturated samples, so the applied effective pressure values can be considered as overburden stress. During the measurements effective pressure was increased from lower to higher values (up-going cycles). The determination of wave arrivals and overall evaluation of results was done using specialized software [26]. of the sample. The correction parameters correspond to the travel time of waves through copper plugs with constant P-and S-wave propagation velocities. The research was performed on 8 samples of core plugs with a diameter of 1″ (2.54 cm) from two boreholes which penetrated the Bjelovar Depression limestone gas reservoir ( Figure 3). The depth range of the investigated core samples is from 805.24 to 974.11 m below the surface which is taken into account during the definition of the applied values of effective pressure.
Prior to the measurement the samples were subjected to a cleaning process which involved organic solvent extraction and drying. Following the sample preparation [25], the effective porosity of the plugs was determined with the Helium expansion method at different confining stress values, namely 5, 7, 9, 11, and 13 MPa.
The ultrasonic velocity measurements were performed on the dry samples in the first place, then the plugs were vacuum saturated with 10 g/L NaCl solution which corresponds to the reservoir brine salinity. The tests were performed at multiple isostatic confining pressure values in both cases, and the applied stresses were identical with the ones used in the porosity determination. A pore pressure was not applied during measurements of the saturated samples, so the applied effective pressure values can be considered as overburden stress. During the measurements effective pressure was increased from lower to higher values (up-going cycles). The determination of wave arrivals and overall evaluation of results was done using specialized software [26].   of the sample. The correction parameters correspond to the travel time of waves through copper plugs with constant P-and S-wave propagation velocities. The research was performed on 8 samples of core plugs with a diameter of 1″ (2.54 cm) from two boreholes which penetrated the Bjelovar Depression limestone gas reservoir ( Figure 3). The depth range of the investigated core samples is from 805.24 to 974.11 m below the surface which is taken into account during the definition of the applied values of effective pressure.
Prior to the measurement the samples were subjected to a cleaning process which involved organic solvent extraction and drying. Following the sample preparation [25], the effective porosity of the plugs was determined with the Helium expansion method at different confining stress values, namely 5, 7, 9, 11, and 13 MPa.
The ultrasonic velocity measurements were performed on the dry samples in the first place, then the plugs were vacuum saturated with 10 g/L NaCl solution which corresponds to the reservoir brine salinity. The tests were performed at multiple isostatic confining pressure values in both cases, and the applied stresses were identical with the ones used in the porosity determination. A pore pressure was not applied during measurements of the saturated samples, so the applied effective pressure values can be considered as overburden stress. During the measurements effective pressure was increased from lower to higher values (up-going cycles). The determination of wave arrivals and overall evaluation of results was done using specialized software [26].   During the implementation of the program, the arrival time of the P-and S-wave was determined at each pressure step ( Figure 4). Based on the sample dimensions, mass, and density, the device automatically calculates the P-and S-wave velocity, Vp/Vs ratio, and other dynamic elastic parameters such as Poisson's ratio, Young's modulus, shear modulus, Lamé constant and acoustic impedance. During the implementation of the program, the arrival time of the P-and S-wave was determined at each pressure step ( Figure 4). Based on the sample dimensions, mass, and density, the device automatically calculates the P-and S-wave velocity, Vp/Vs ratio, and other dynamic elastic parameters such as Poisson's ratio, Young's modulus, shear modulus, Lamé constant and acoustic impedance.

Results and Interpretation
The P-and S-wave velocities of the Bjelovar Depression limestones were measured as a function of effective pressure under dry and saturated conditions ( Figure 5). An analysis of the measurement results (Table 1) was made with respect to saturation, porosity and effective stress. The correlation of P-wave and S-wave velocities and elastic parameters were also analyzed.

Results and Interpretation
The P-and S-wave velocities of the Bjelovar Depression limestones were measured as a function of effective pressure under dry and saturated conditions ( Figure 5). An analysis of the measurement results (Table 1) was made with respect to saturation, porosity and effective stress. The correlation of P-wave and S-wave velocities and elastic parameters were also analyzed. During the implementation of the program, the arrival time of the P-and S-wave was determined at each pressure step ( Figure 4). Based on the sample dimensions, mass, and density, the device automatically calculates the P-and S-wave velocity, Vp/Vs ratio, and other dynamic elastic parameters such as Poisson's ratio, Young's modulus, shear modulus, Lamé constant and acoustic impedance.

Results and Interpretation
The P-and S-wave velocities of the Bjelovar Depression limestones were measured as a function of effective pressure under dry and saturated conditions ( Figure 5). An analysis of the measurement results (Table 1) was made with respect to saturation, porosity and effective stress. The correlation of P-wave and S-wave velocities and elastic parameters were also analyzed.

Correlation of Velocity and Porosity
The measured velocities on dry and water saturated samples are correlated to porosity and presented in Figure 6. The porosity under different load pressures is the greatest Laboratory results were also interpreted and verified using two theoretical rock physics models, the Kuster and Toksöz model [27] and the Xu-Payne model [28]. The Kuster and Toksöz (1974) model presented randomly oriented pores with the assumption of low inclusion concentration. The model implies only one pore type and in the present study results are interpreted assuming penny shaped pores. Pore shape is defined by pore aspect ratio (α) which represents ratio of the perpendicular axis of the pore (Figures 7 and  8).
Xu and Payne, (2009) presented a rock physics model for carbonates based on the Xu White model [29] for clastic rocks. The rock physics model includes three pore types, characteristic for carbonate rocks: moldic pores (the rounded stiff pores), interparticle pores and microcracks.
Xu and Payne, (2009) [28] assumed that the interparticle pore system is the most common for carbonates and pore aspect ratio trend of 0.15 presents reference velocityporosity trendline. Lines below the reference trendline present the pore system with the increasing volume of microcracks with best fit aspect ratio of 0.02 while lines above the reference trendline present the increasing volume of stiff pores with an aspect ratio of 0.8 ( Figure 8). Although the velocities decrease uniformly with increasing porosity regardless of the effective stress, the velocities oscillate with increasing effective stress. The reason for these oscillations is the sudden decrease in porosity during the application of the lowest effective stress and the activation of the predefined crack system with a further increase in the effective stress, e.g., some samples at 13 MPa have lower velocities and higher porosity values than expected (Figure 6a,b).
Laboratory results were also interpreted and verified using two theoretical rock physics models, the Kuster and Toksöz model [27] and the Xu-Payne model [28]. The Kuster and Toksöz (1974) model presented randomly oriented pores with the assumption of low inclusion concentration. The model implies only one pore type and in the present study results are interpreted assuming penny shaped pores. Pore shape is defined by pore aspect ratio (α) which represents ratio of the perpendicular axis of the pore (Figures 7 and 8). trendline present the increasing volume of stiff pores with an aspect ratio of 0.8 ( Figure  8).  The Kuster and Toksöz model indicates penny shaped pore distribution between aspect ratios 0.04 and 1 (Figure 7). The results are mainly constrained between 0.06 and 0.08 values of aspect ratio. Figure 7 shows that the Kuster and Toksöz model cannot describe the data well enough since it is not possible to use a unique pore aspect ratio for all results.
On the other hand, the interpretation of the study results with the Xu-Payne model enables the definition of the pore system of the studied Pannonian Basin limestone. In  trendline present the increasing volume of stiff pores with an aspect ratio of 0.8 ( Figure  8).  The Kuster and Toksöz model indicates penny shaped pore distribution between aspect ratios 0.04 and 1 (Figure 7). The results are mainly constrained between 0.06 and 0.08 values of aspect ratio. Figure 7 shows that the Kuster and Toksöz model cannot describe the data well enough since it is not possible to use a unique pore aspect ratio for all results.
On the other hand, the interpretation of the study results with the Xu-Payne model enables the definition of the pore system of the studied Pannonian Basin limestone. In Xu and Payne, (2009) presented a rock physics model for carbonates based on the Xu White model [29] for clastic rocks. The rock physics model includes three pore types, characteristic for carbonate rocks: moldic pores (the rounded stiff pores), interparticle pores and microcracks.
Xu and Payne, (2009) [28] assumed that the interparticle pore system is the most common for carbonates and pore aspect ratio trend of 0.15 presents reference velocityporosity trendline. Lines below the reference trendline present the pore system with the increasing volume of microcracks with best fit aspect ratio of 0.02 while lines above the reference trendline present the increasing volume of stiff pores with an aspect ratio of 0.8 ( Figure 8). The Kuster and Toksöz model indicates penny shaped pore distribution between aspect ratios 0.04 and 1 (Figure 7). The results are mainly constrained between 0.06 and 0.08 values of aspect ratio. Figure 7 shows that the Kuster and Toksöz model cannot describe the data well enough since it is not possible to use a unique pore aspect ratio for all results.
On the other hand, the interpretation of the study results with the Xu-Payne model enables the definition of the pore system of the studied Pannonian Basin limestone. In Figure 8, the P-wave velocity is plotted versus porosity for both dry and saturated samples. The crossplot indicates that the interparticle pore system dominates with a fraction of microcracks since all the data fall below the reference line. Most of the data points are grouped below the 20% crack pore volume line, so according to the interpretation of the results in the Xu-Payne model the fraction of microcracks is between 20% and 35% of the total pore space.
The interpretation is in accordance with the velocity and porosity results and supports the explanation of the presence of predominant fractures as the cause of the velocity oscillation with increasing effective stress.

Correlation of P-Wave and S-Wave Velocities
Measured P-and S-wave velocities were also compared with each other (Figure 9). Generally, the correlation of P and S velocities is significant for quantitative interpretation and characterization of reservoirs and many authors have published empirical relations for different lithologies, e.g., [30].
Energies 2021, 14, x FOR PEER REVIEW 9 of 18 Figure 8, the P-wave velocity is plotted versus porosity for both dry and saturated samples. The crossplot indicates that the interparticle pore system dominates with a fraction of microcracks since all the data fall below the reference line. Most of the data points are grouped below the 20% crack pore volume line, so according to the interpretation of the results in the Xu-Payne model the fraction of microcracks is between 20% and 35% of the total pore space. The interpretation is in accordance with the velocity and porosity results and supports the explanation of the presence of predominant fractures as the cause of the velocity oscillation with increasing effective stress.

Correlation of P-Wave and S-Wave Velocities
Measured P-and S-wave velocities were also compared with each other (Figure 9). Generally, the correlation of P and S velocities is significant for quantitative interpretation and characterization of reservoirs and many authors have published empirical relations for different lithologies, e.g., [30].
The results indicate that P-wave velocities of dry samples are lower than velocities of saturated ones, from 15 m/s to 330 m/s or 0.3% to 7.7% (Figure 9a). S-wave velocities of dry samples are lower than saturated (Figure 9b). Although theoretically S-wave velocities should not change during fluid replacement, changes in velocity values are due to density change of the medium by replacing the pore fluid, air to water with 10 g/L NaCl, and moreover because of shear modulus weakening or strengthening. The results show that S-wave velocities of saturated limestone samples are 3 m/s to 290 m/s or 0.1% to 11.6% lower than dry samples. Since the velocity-porosity crossplot indicates uniform velocity trend regardless of effective pressure, correlation between Vp and Vs is analyzed together for all values of effective pressure (Figure 10). The results indicate that P-wave velocities of dry samples are lower than velocities of saturated ones, from 15 m/s to 330 m/s or 0.3% to 7.7% (Figure 9a). S-wave velocities of dry samples are lower than saturated (Figure 9b). Although theoretically S-wave velocities should not change during fluid replacement, changes in velocity values are due to density change of the medium by replacing the pore fluid, air to water with 10 g/L NaCl, and moreover because of shear modulus weakening or strengthening. The results show that S-wave velocities of saturated limestone samples are 3 m/s to 290 m/s or 0.1% to 11.6% lower than dry samples.
Since the velocity-porosity crossplot indicates uniform velocity trend regardless of effective pressure, correlation between Vp and Vs is analyzed together for all values of effective pressure (Figure 10).

Elastic Parameters
Besides seismic velocities elastic moduli were determined for all dry and water saturated samples ( Table 2). All analyzed samples show an increase in bulk modulus after water saturation which is expected according to many earlier published results, e.g., [5,6,10]. The maximum increase in value is 32.9% and the minimum is 8.6%. The results indicate that the bulk modulus increase with increasing effective stress (Figure 11a).

Elastic Parameters
Besides seismic velocities elastic moduli were determined for all dry and water saturated samples ( Table 2). All analyzed samples show an increase in bulk modulus after water saturation which is expected according to many earlier published results, e.g., [5,6,10]. The maximum increase in value is 32.9% and the minimum is 8.6%. The results indicate that the bulk modulus increase with increasing effective stress (Figure 11a).  On the other hand, the effect of shear modulus weakening with water saturation is visible on the studied samples (Figure 11b). The greatest shear modulus weakening of the saturated samples is 19% for sample No. 2. But some samples such as No. 4, No. 7 and No. 8 show an almost linear (theoretical trend). These samples have a slight increase in the shear modulus of saturated samples from 0.20 to 1.3%. Although, their porosity is low with values from 8 to 10%, the relationship between porosity and shear modulus increase was not observed because other samples with the same porosity show the effect of shear modulus weakening. Results in Figure 11b show that three samples have a shear modulus strengthening effect, with pronounced secondary porosity which is consistent with the conclusions of Baechle et al. (2009) [5].

Sample No. CP (MPa) K (GPa) G (GPa) E (GPa) ν (-) K (GPa) G (GPa) E (GPa) ν (-)
The values of the Poisson's ratio vary with porosity and depend on the state of saturation, i.e., water saturated or dry sample (Figure 12). The Poisson ratio of dry samples is between 0.203 and 0.269 (Table 2) generally showing decreasing trend with increasing porosity. For saturated samples the trend is opposite presenting increasing Poisson's ratio with porosity ( Figure 12). The Poisson ratio of water saturated samples is 0.273-0.332 ( Table 2). The variance of the Poisson's ratio with effective stress is not significant because the Poisson's ratio depends on the Vp/Vs values which also do not show dispersion with the change of effective stress. The Vp/Vs values of dry samples decrease slightly, while those of saturated samples increase slightly with increasing porosity. The mean value of Vp/Vs of dry samples is 1.72 and saturated 1.86 (Table 1).  (Table 1).
(a) (b) Figure 11. Crossplot of elastic moduli between dry and saturated samples: (a) bulk modulus and (b) shear modulus. Figure 11. Crossplot of elastic moduli between dry and saturated samples: (a) bulk modulus and (b) shear modulus.
To understand and interpret the elastic property of the samples, we performed Gassmann fluid substitution using directly the dry rock results and compared them to results of saturated samples. Gegenhuber (2015) [31] presented the direct implication of dry rock laboratory results into the Gassmann formula for Austrian carbonates.
The bulk modulus resulting from the test on dry samples (Table 1) was used directly in the Gassmann equation. This method assumes uniform pore pressure in the rock [32]. The calculation on dry rock gives an overall good agreement with the saturated samples, but the Gassmann equation gives slightly overestimated results. Figure 13 presents detailed values of the difference between modeled and measured bulk modulus, and the highest value of overestimation is 4 GPa.  To understand and interpret the elastic property of the samples, we performed Gassmann fluid substitution using directly the dry rock results and compared them to results of saturated samples. Gegenhuber (2015) [31] presented the direct implication of dry rock laboratory results into the Gassmann formula for Austrian carbonates. The bulk modulus resulting from the test on dry samples (Table 1) was used directly in the Gassmann equation. This method assumes uniform pore pressure in the rock [32]. The calculation on dry rock gives an overall good agreement with the saturated samples, but the Gassmann equation gives slightly overestimated results. Figure 13 presents detailed values of the difference between modeled and measured bulk modulus, and the highest value of overestimation is 4 GPa.

Discussion
Current study results help to understand the elastic properties of limestones in the southern Pannonian Basin, based on core samples from the Bjelovar Depression, which is necessary for any process of reservoir characterization (e.g., seismic, geomechanics).
The results of the Vp−Vs velocity relationship can be compared with published papers and results. So far, the most detailed analysis of carbonate ultrasonic velocities has been published by Rafavich et al. (1984) [2] based on the laboratory analysis of cores from 4 wells and a total of 93 samples, but here only limestone samples were considered. Ultrasonic velocities were analyzed under an effective stress of 6000 psi, or 40 MPa, on water-saturated samples. Study of Assefa et al. (2003) [33] was performed on limestone core samples from Southern England under in situ conditions at confining pressure of 50 MPa and demonstrated generally slightly lower both P-and S-wave velocities. Figure 15 shows a comparison of Vp and Vs velocities and porosity based on the results of Rafavich et al. (1984) [2], Assefa et al., (2003) [33] with the results of current research. Rafavich et al. (1984) [2] concluded that the major influence on velocity has porosity and bulk density, while the influence of calcite and dolomite is minor compared to the first two parameters. Our observations show the same trend of decreasing velocities with increasing porosity. It is important to emphasize that there is a difference in laboratory measurements conditions. The results of Rafavich et al. (1984) [2] are based on cores from a depth range of 1200 to 3100 m from the surface with confining pressure of 40 MPa, while the results of the current study are based on cores from a depth range of 805.24 to 974.11 m from the surface with different effective pressures.
ity and bulk density, while the influence of calcite and dolomite is minor compared to the first two parameters. Our observations show the same trend of decreasing velocities with increasing porosity. It is important to emphasize that there is a difference in laboratory measurements conditions. The results of Rafavich et al. (1984) [2] are based on cores from a depth range of 1200 to 3100 m from the surface with confining pressure of 40 MPa, while the results of the current study are based on cores from a depth range of 805.24 to 974.11 m from the surface with different effective pressures. Such consistency in velocity-porosity trends have been discussed by [34] on studied wireline data of chalk despite very large distance to reservoir location. They conclude that diagenetic trends are universal in carbonates. If we consider also results given by [5,7,33] correlation of Vp and Vs ( Figure 16) with different confining pressure have credibility for defining the empirical relation defined with formula 1 and 2. Such consistency in velocity-porosity trends have been discussed by [34] on studied wireline data of chalk despite very large distance to reservoir location. They conclude that diagenetic trends are universal in carbonates. If we consider also results given by [5,7,33] correlation of Vp and Vs ( Figure 16) with different confining pressure have credibility for defining the empirical relation defined with formula 1 and 2.   [5] imply that Gassmann theory [8] does not consider a constant shear modulus with fluid substitution, but this is considered by numerous Gassmann assumptions for porous media. Therefore, using the Figure 16. Comparison of Vp-Vs correlation results in the current study with previous results of [2,5,7,33] obtained for different effective pressures. Baechle et al., (2009) [5] have studied elastic properties of limestone samples under dry and water saturated conditions with confining pressure from 2 to 80 MPa, but they presented results only for 40 MPa. They concluded that velocities have an inverse relationship to porosity. Elastic properties and prediction models for fully saturated carbonate from Saudi Arabia and Canada are presented by Bakhorji, (2010) [7] study. Laboratory measurements were obtained with different effective pressure steps (2-25 MPa) which overlaps with effective pressure of the current study.
The results of Baechle et al. (2009) [5] and Assefa et al. (2003) [33] indicated decrease in shear modulus for saturated samples. Baechle et al. (2009) [5] imply that Gassmann theory [8] does not consider a constant shear modulus with fluid substitution, but this is considered by numerous Gassmann assumptions for porous media. Therefore, using the Gassmann equation to estimate velocities in these samples may be inaccurate. Zhao et al. (2021) [35] have implemented an extended Gassmann equation for heterogeneous rocks with complex properties. They concluded that microcracks and patchy saturation cause attenuation. Regnet et al. (2015) [32] state that the use of dispersion analysis excludes external source, such as cracks due to mechanical compaction and emphasizes the physical relationship. Zhao et al. (2013) [36] assume a random orientation of cracks which significantly affects the seismic wave, especially in the presence of cracks with low aspect ratio. Also, if rocks have two pore systems, cracks and stiff pores, velocity dispersion may occur as a result of elastic pore heterogeneity.
Dispersion affects the velocity-porosity relationship, but according to our results the velocity-porosity relationship is uniform, even for samples with higher dispersion values (2 and 4). Still, considering the defined two-pore system and the presence of stylolite filled with clay minerals, we cannot exclude the influence of anisotropy on our results.
Although, this research gives insight into limestone properties, it also gives an opportunity for further research. Elastic parameters presented in this study can be correlated with static elastic parameters performed on the same samples. Such results would provide correction factors that can be applied on well logs directly and would acquire scientific and industrial benefits (e.g., reservoir workovers as fracking).

Conclusions
In the lack of published limestone elastic properties in the Pannonian Basin, especially in Croatia, this paper presents results of ultrasonic velocity measurements, porosity measurements and dynamic elastic properties determination carried out for the first time on the limestone core samples in this part of the Pannonian Basin. The research was performed on 8 core plugs from two boreholes in the Bjelovar Depression, which were in two states, dry and brine saturated. The tests were done at multiple isostatic confining pressure values in both cases.
The porosity under different load pressures range from 5.79% to 15.06%. The variation of velocities with respect to porosity shows the decreasing trend of P-and S-wave velocities with increasing porosity regardless of the effective stress, in a dry and saturated states. Although the velocities decrease uniformly with increasing porosity some samples show velocity scattering with increasing effective stress. The reason for these oscillations is a sudden decrease in porosity during the application of the lowest effective stress and the activation of a predefined crack system with further increase in the effective stress. Based on the Xu-Payne model, the pore type of limestone samples is defined as interparticle pore system with a fraction of cracks from 20% to 35%.
The results show that the P-wave velocities of dry samples are 0.3% to 7.7% lower than those of saturated samples, while the S-wave velocities of saturated limestone samples are 0.1% to 11.6% lower than those of dry samples. The results indicate uniform trend of Vp and Vs for all values of confining pressure, and are consistent with the results in similar published papers [2,5,7,33]. Therefore, an empirical relationship between Vp and Vs for dry and brine saturated limestones has been defined.
All analyzed samples show an increase in bulk modulus after water saturation with a maximum increase in value of 32.9% and a minimum increase of 8.6%. Though, the effect of shear modulus weakening due to water saturation is observed. The greatest shear modulus weakening of the saturated sample is 19% for sample No. 2, although some samples show an almost linear (theoretical) trend.
We performed Gassmann fluid substitution using directly the dry rock results and compared them to the results of saturated samples. The calculation results based on dry rock gives an overall good fit to the saturated samples, but the Gassmann equation gives slightly overestimated results.
The major benefits of this research are the first insights on elastic properties of the Pannonian Basin limestones which provides the basis for further study beneficial for science and industry. The defined empirical relationship allows transformation of Vp to Vs from well log data considering the modified Gassmann workflow.