Temperature-Dependent Wave Velocities of Heavy Oil-Saturated Rocks

Understanding the dependence of the rock properties on temperature is essential when dealing with heavy oil reservoirs. Reported rock physics models can hardly capture the effect of temperature on wave velocities. We propose a dual-porosity temperaturedependent model based on the coherent potential approximation, combining temperatureand frequency-dependent empirical equations for pore fluids with the David-Zimmerman model. The Maxwell model is adopted to obtain the complex shear modulus as a function of temperature and frequency. To verify the validity of the model, a glycerol-saturated tight sandstone and three heavy oil sand samples are considered. The comparison between the predicted and measured wave velocities shows that the model can quantitatively describe the behavior as a function of temperature. We find that there is a viscosity threshold (liquid point), where the Pand S-wave velocity variation trends change, while the porosity has no effect.


Introduction
Reservoir rocks can be considered as heterogeneous porous viscoelastic media, fully or partially saturated with fluids which have a dissimilar behavior as a function of depth, temperature, and pressure. Seismic wave attenuation and velocity dispersion are sensitive to the fluid types and corresponding environmental conditions [1,2] and can be used to obtain information about the rock properties. Then, it is important to investigate the influence of the saturating fluid on these properties and the fluid-skeleton interaction.
Fluids have very dissimilar properties. For example, the thermal properties of light gases with low density, modulus, and velocity are quite distinct from those of heavy oil with a high viscosity, as well as those of brine and light oil. Fluid properties can be estimated as a function of the local pressure and temperature conditions, using empirical relations [3][4][5][6]. A large number of theoretical and experimental studies have been performed on the elastic wave velocities in rocks under different environmental (temperature and pressure) conditions [7][8][9][10].
However, modeling on the physical properties of heavy oil is difficult due to its high viscosity. Han et al. [11] pointed out that the USGS (United States Geological Survey) defines heavy oil as a dense viscous oil chemically characterized by the asphaltene content. The API (American Petroleum Institute of Oil Gravity) of heavy oil is defined as 22 to 10 or less (ultraheavy oil or bitumen). Many authors have studied the elastic properties of heavy oil-saturated rocks by means of laboratory measurements [12][13][14][15]. Heavy oil is solid at room temperature and pressure but melts when heated, thus decreasing the shear wave velocity [16]. Appropriate viscoelastic models can be adopted to obtain the properties of heavy oil-saturated rocks [11,[17][18][19][20][21][22], and the shear modulus of the fluid cannot be neglected in this case [23]. Han et al. [20] show that heavy oil becomes relatively dispersed when the temperature is below 60°C and the predicted values obtained with the Gassmann-Biot theory do not match the measured experimental data [24,25].
The conventional fluid substitution method is not applicable here due to the high viscosity of heavy oil. The problem has been investigated by using the CPA (coherent potential approximation) method, which takes into account the effect of viscosity and temperature and considers the fluid as a viscoelastic medium [18,21,26]. However, there are still uncertainties in the proposed models, which limit the practical applications.
We develop a temperature-dependent rock physics model based on the double-porosity CPA method and the micropore structure theory proposed by David & Zimmerman [27]. The CPA method allows us to consider a pore fluid of high viscosity described by the Maxwell viscoelastic model. The Maxwell model has been successfully used to describe molten material in the Earth by Carcione et al. [28]. Then, the temperature-dependent wave velocities are compared to experimental data.

A Temperature-Dependent Rock
Physics Model 2.1. Single-Porosity Coherent Potential Approximation Medium (CPA) Method. The conventional modeling approach of fluid substitution is not suitable for heavy oilsaturated rocks, because the fluid properties render the main assumptions of the Gassmann-Biot theory invalid [18]. Heavy oil exhibits Newtonian fluid properties at high temperatures and resembles elastic solids at room temperature. In order to analyze the influence of high-viscosity fluids, we consider the CPA method. When the fluid saturation is low, the solid-fluid mixture is solid with a specific shape of the fluid inclusions. On the other hand, when the solid content is low, the mixture is a suspension of solid grains in the fluid [18,21,26]. The equations are where ϕ is the porosity, K and G are the bulk and shear moduli of the rock, respectively, K f and G f are the bulk and shear moduli of the pore fluid, respectively, K s and G s are the bulk and shear moduli of the rock matrix, respectively, and P and Q are shape factors [26,29].

Fluid Model.
Heavy oil is a viscoelastic medium with three phases (fluid, quasisolid, and solid) and two temperature points (liquid point and glass point) [11]. When the temperature is below the glass point, the fluid is considered part of the solid frame due to its high shear viscosity. At high temperatures, the viscosity is low and can be treated as a conventional fluid without effects on the rock frame. Between these two phases, the fluid presents a viscoelastic state, where the waves exhibit strong dispersion and attenuation. It is in the quasisolid phase with a finite shear modulus. By using the CPA, in addition to considering the properties of the fluid, i.e., shear viscosity ηðTÞ, density ρ f ðTÞ, and bulk modulus K f ðTÞ, the shear modulus of the fluid G f ðTÞ should be incorporated. Since the properties of viscoelastic materials are frequency dependent, it is reasonable to consider them in the frequency domain [30][31][32]. Here, the Maxwell model is adopted to obtain the complex shear modulus as a function of temperature and frequency [33]: where G ∞ is the modulus at high frequencies, τ is a relaxation time, ω is the frequency, i = ffiffiffiffiffi ffi −1 p , and ηðTÞ is the viscosity of fluid varying with temperature.
By substituting equation (3) into (2), we obtain where Re ðG f ðTÞÞ and lmðG f ðTÞÞ are the storage and loss moduli, respectively [34]. The fluid shear modulus considered in this study refers to the storage modulus.
2.3. Microscopic Pore Structure Model. Makarynska et al. [21] showed that the use of the single-porosity CPA method is insufficient to explain the effects of temperature and frequency on heavy oil-saturated rocks. The reason is that the stiff pores occupy most of the pore but there are also remaining cracks or grain contacts, which affect the elastic moduli of the rock as a function of temperature and pressure. We consider the double-porosity CPA method developed by Makarynska et al. [21], but here, we also consider the effects of pressure and temperature on the cracks, absent in the previous approach.
The crack porosity is [27,35] where β is the crack density and α p is the pore aspect ratio. Firstly, ultrasonic experimental data varying with effective pressure at different temperatures are required. We consider that the cracks tend to close at a high effective pressure and only stiff pores are present, whose moduli can be estimated based on the differential effective medium theory [36]. According to this scheme, the effective bulk and shear moduli as a function of the crack density are [37] determined.
where v stiff = ð3K stiff ,dem − 2G stiff ,dem Þ/ð6K stiff ,dem + 2G stiff ,dem Þ is the Poisson ratio of the stiff pores and K stiff ,dem and G stiff ,dem are the related moduli. Then, the crack density at a 2 Lithosphere given temperature and effective pressure can be obtained by least-square regression based on equations (7)- (9). The quantitative relationship between the pore aspect ratio α p and the effective pressure p is established: where E s is the Young modulus at high effective pressures. It is defined as E s = 3K s ½1 − 2ν s , where ν s is the Poisson ratio. Substituting equations (7)-(10) into (6), the crack porosity varying with pressure at the same temperature can be obtained. Next, the discrete data points of crack porosity at different temperatures and the same effective pressure must be obtained. This relation between crack porosity and temperature is then substituted into equation (1).

Temperature-Dependent Double-Porosity CPA. Equation
(1) can be simplified. In the CPA method with m inclusions, the equations are [28,38]: Using this method, we consider the inclusion and background, and in addition, the pore space consists of stiff pores and cracks. Thus, we have three components, i.e., the solid material, the stiff porosity that occupies most of the pore space, and the crack porosity at grain contacts.
By substituting equations (2)-(10) into equation (11), the temperature-dependent CPA model with two pore phases is obtained as where ϕ 1 is the stiff porosity and ϕ = ϕ 1 + ϕ 2 . Coupled with the implicit equation (12) for K and G, the effective bulk modulus K n+1 ðTÞ and shear modulus G n+1 ðTÞ of the saturated rock can be computed by numerical iteration according to the following scheme as The iterative process requires the initial K 1 and G 1 , which can be computed by using the Voigt-Reuss-Hill (V-R-H) average [39][40][41][42]. P n and Q n are evaluated in this process by using the nth approximations to the moduli of rock. We substitute the initial values into equation (12) and iterate n + 1 times, thus obtaining K n+1 ðTÞ and G n+1 ðTÞ.
Then, the P-and S-wave velocities and density of the saturated rock are where ρ s is the matrix density.
Here, we consider five temperature-dependent parameters, namely, viscosity ηðTÞ, density ρ f ðTÞ, bulk modulus K f ðTÞ, shear modulus G f ðTÞ, and crack porosity ϕ 2 ðTÞ. When the temperature is low, the effect of the bulk viscosity on the bulk modulus is less than that of the shear viscosity on the shear modulus [18], the bulk viscosity can be neglected [43,44]. As the temperature increases, the shear modulus can be neglected. Figure 1 shows a diagram of the temperature-dependent model, and Figure 2 gives the detailed workflow. We first estimate the bulk modulus K 1 and shear modulus G 1 of the rock matrix based on the mineral content. Then, the crack porosity ϕ 2 ðTÞ is obtained by equations (6)- (10). The thermophysical properties of the fluids are given by the Batzle-Wang empirical formulae [3]. The storage shear modulus Re ðG f ðTÞÞ is obtained by equations (4) and (5). Next, these parameters are substituted into equation (1) of the single-porosity CPA model to obtain equation (12). Finally, equation (13) is solved by iterations, where K 1 and G 1 are assumed to be the initial bulk and shear moduli of the rock matrix, respectively. Finally, equation (14) is used to obtain the P-and S-wave velocities of rock saturated with heavy oil.

Example: Fluid Properties
We consider low-and high-viscosity fluids, namely, water and Uvalde heavy oil (API = −5°), respectively [7]. We take a pure quartz sandstone as an example, whose properties are listed in Table 1. The porosity is 35%.
Based on the temperature-dependent double-porosity CPA, we analyze water-saturated and heavy oil-saturated sandstones separately. The bulk modulus of Uvalde heavy oil is 3 GPa [45], and the corresponding shear modulus is obtained from equation (4). The bulk modulus of water can be obtained according to Batzle & Wang [3] at the same temperature. The porosity is 35% and the frequency is 1 MHz. Then, we use equation (12) to obtain the moduli of the saturated rock. Figure 3 shows the bulk and shear moduli as a function of porosity and different stiff pore aspect ratios. At the same aspect ratio, the moduli of the heavy oilsaturated sandstone are higher than those of the watersaturated sandstone. With an increasing aspect ratio, the moduli also increase. Figure 4 shows the same moduli as a function of temperature. As the temperature increases, the bulk and shear moduli of Uvalde heavy oil significantly decrease. In contrast, when the temperature exceeds 50°C, the bulk modulus decreases in a lesser extent, whereas the shear modulus remains almost constant. This transition behavior in the lower temperature range is associated with the fact that the viscosity affects the bulk and shear moduli of heavy oil [7,44]. Since the shear modulus of water is zero, the shear modulus of the water-saturated rock generally does not depend on temperature. On the other hand, Uvalde heavy oil is non-Newtonian fluid with high viscosity, the quasisolid phase at low temperatures, and its shear modulus is not neg-ligible. It gradually approaches a fluid phase as temperature increases, where the shear modulus vanishes. The modulus of the heavy oil-saturated sandstone with temperature is nearly linear in the high-temperature range, and the aspect ratio of the stiff pores has a noticeable effect.
In addition to the effect of temperature, the frequency also affects the modulus of the fluid. Next, we consider an aspect ratio of the stiff pores of 0.5 and the properties of Table 1. Equation (12) gives the moduli of the saturated rock. Figure 5 shows the moduli as a function of frequency and temperature. At 40°C, there is dispersion, and when the temperature increases from 40°C to 60°C, the moduli dispersion sharply decreases. At the same temperature, the   Lithosphere moduli of Uvalde heavy oil significantly increase with frequency. When the temperature is such that the infill is in a quasisolid state, viscosity has an important effect and causes modulus dispersion due to its frequency dependence.

Results and Discussions
4.1. Examples. We consider ultrasonic experimental data reported in the literature, and for each sample, at least three temperature ranges are included. Table 2 lists the rock and  fluid types, while Tables 3 and 4 show the properties.
We plot the normalized velocities as a function of temperature (V P /V P0 and V S /V S0 , where V P0 and V S0 are the P-and S-wave velocities at room temperature, respectively). Figure 6(a) shows the results for the glycerol-saturated tight sandstones reported by Yin et al. [46] at different temperatures and an effective pressure of 7 MPa. The velocities decrease with temperature, as expected. Figure 6(b) shows the results for the three oil sand samples with high porosity by Li et al. [44] at different temperatures with a confining pressure of 1100 psi. We can see that the decreasing trend of the S-wave velocity is more pronounced than that the of P-wave velocity.

Laboratory Measurements in Glycerol-Saturated Tight
Sandstone. By considering sample T1 whose properties are given in Table 3, we used the proposed model to analyze the relation between the temperature and the P-and Swave velocities. The tight sandstone sample was collected from the northeastern Sichuan Basin, China [47]. The main mineral components include quartz, feldspar, and calcite, along with many other interstitials such as clay and flaky silicate. The sample has a diameter of 38 mm and a length of 70 mm. For the glycerol-saturated condition, ultrasonic experiments were carried out at 23°C, 40°C, and 60°C, with a frequency of 1 MHz. During the experiment, the pore pressure was set at 1 MPa. The specific properties are shown in Table 3.

Effect of Temperature on the Properties of Glycerol.
Since glycerol is a high-viscosity fluid, its thermal properties change with temperature [48] as shown in Figure 7, especially when the temperature is close to the liquid and glass points. Figure 7(a) shows that the viscosity of glycerol decreases with increasing temperature. The liquid-point temperature can be defined when the viscosity of the fluid rapidly drops to 10 3 cP. Below the liquid point (25.5°C), the glycerol is a quasisolid phase. As the temperature exceeds the liquid point, glycerol is a fluid. Figure 7(b) shows the density and bulk modulus of glycerol as a function of temperature. At a given pressure, the bulk modulus of glycerol decreases almost linearly with increasing temperature. The lower the temperature, the higher the bulk modulus, since glycerol becomes stiffer. In addition, the density of glycerol decreases rapidly with increasing temperature. Figure 8 shows the shear modulus of glycerol computed with the viscoelastic Maxwell model (equations (4) and (5)). Each curve corresponds to different frequencies. The infinite shear modulus (G ∞ ) decreases linearly with temperature, which is obtained by a linear fitting [49]. At ultrasonic frequencies, the shear modulus of glycerol decreases to zero approximately at the liquid point. As the frequency increases, the liquid point moves to the right. The velocity dependence on frequency occurs when a high-viscosity fluid is in the quasisolid phase. In fact, the high-viscosity fluid of the glass phase is related to the high viscosity (>10 15 cP), which is generally considered as elastic, whereas the molecules in this medium are stable and equivalent to a rigid material, similar to crystalline solids. High-viscosity fluids in the fluid phase can also be considered elastic. Therefore, the effect of frequency on velocity is coupled with viscosity [7]. The temperature and viscosity of the liquid point depend on frequency. It can be considered as a relaxation phenomenon, i.e., the effective stiffness depends on the rate of shear modulus deformation of the pore fluid. It is shown in Figure 7(a) that at 1 MHz, the viscosity at the glycerol liquid point is about 1000 cP. However, at 1000 Hz, the temperature of the glycerol liquid point is lower than that at 1 MHz, so that the viscosity at the liquid point is higher. Therefore, for a high-viscosity fluid-saturated rock, it is very important to consider the influence of temperature and frequency on the shear viscosity of the fluid.

The
Effect of Temperature on Velocity. Following the workflow of Figure 2, we obtain the crack porosity according to the V P and V S provided by Yin et al. [46]. Figure 9 shows the crack porosity as a function of temperature. It decreases with effective pressure at a given temperature. At low temperatures, glycerol is nearly solid and cannot flow in the cracks and the grain contact is stiff. However, with the increase of temperature, the viscosity of glycerol decreases     7 Lithosphere rapidly and the fluidity increases. Therefore, in addition to the thermal properties of the minerals, the thermal expansion of the fluid has an important effect on cracks [50,51]. Since the normal stiffness of grain contact depends on the storage modulus of the pore-filling medium, the cracks may also lead to the dispersion of P-and S-waves in rocks saturated with a high-viscosity fluid [52]. Figure 10 shows the results predicted by the proposed model, which are in good agreement with the measured Pand S-wave velocities from ultrasonic measurements, indicating that the model can effectively describe the variation of velocities with temperature in glycerol-saturated tight sandstones. The high-viscosity fluid highly affects the properties of the saturated rock. When the temperature is less than the glass point, the fluid can be considered as part of the solid frame. When the temperature is higher than the liquid point, the fluid can be treated as a conventional fluid (e.g., brine or light oil). Both the P-and S-wave velocities decrease with increasing temperature. The velocities significantly decrease with increasing temperature. As the temperature approaches the liquid point, glycerol is in the fluid phase and the velocities linearly decrease with temperature slightly. In this case, the shear modulus of glycerol is negligible.

Liquid Point.
According to the specific properties of high-viscosity fluids, Han et al. [11] proposed the threshold value of fluid viscosity, namely, the liquid point, to describe the nonlinear P-and S-wave velocity behavior. It is assumed that the normalized temperature dictates the behavior of the high-viscosity fluid,   Lithosphere where T w is the temperature at which the fluid viscosity is equal to the water viscosity (1 cP) and T g refers to the glass point, which is the temperature at 10 15 cP. Figure 11 shows the P-and S-wave velocities of glycerolsaturated tight sandstone as a function of the normalized temperature. If T n is less than 0.58, there is a linear relationship between the P-wave (S-wave) velocity and the normalized temperature. On the other hand, when T n is higher than 0.58, the behavior is nonlinear. Therefore, 0.58 corresponds to the liquid point. For the same fluid, the liquid point (the viscosity threshold) is different at different frequencies (see Figure 8).

Laboratory Measurements in Heavy Oil
Sand. The experimental data of three oil sand samples reported by Li et al. [44] are analyzed. The main mineral components of these samples include quartz and clay. The pore fluid is a heavy oil with an API density of 6.6°. Ultrasonic experiments were carried out at four temperature points with a frequency of 1 MHz. During the experiment, the pore pressure was set to 0 psi. The rock physical properties are given in Table 4.

Effect of Temperature on the Properties of Heavy Oil.
The relation between viscosity and temperature is shown in Figure 12, where we can see that viscosity decreases with     9 Lithosphere temperature. When the viscosity reaches 10 3 cP, it can be considered as the liquid point. When the temperature is below the liquid point (48.7°C), the red curve represents the quasisolid phase. As temperature increases, it exceeds the liquid point and the blue curve represents the fluid phase. Figure 13 shows the relation between the shear modulus and temperature of the heavy oil at different frequencies, according to equations (4) and (5). In this case, the infinite shear modulus (G ∞ ) of heavy oil decreases linearly with temperature and can be obtained by a linear regression fitting of the ultrasonic shear modulus [44]. At ultrasonic frequencies, the shear modulus is zero at approximately 48.7°C. It shows that the liquid point moves to the right side with increasing frequency. At 1 MHz, the viscosity of the heavy oil liquid point is about 1000 cP.

Effect of Temperature on Wave Velocities and Fluid
Point. Figures 14(a)-16(a) show the results obtained with the proposed model, which are in agreement with the velocities measured in ultrasonic experiments of the three oil sands, indicating that the model effectively describes the behavior with temperature. The velocities decrease with increasing temperature, and the trend is the same as that of Figure 6. When the temperature is less than the liquidpoint one, the heavy oil is at a quasisolid phase with a rapidly decreasing viscosity and reacts as a mobile pore fluid with low shear viscosity, so the velocities decrease sharply with increasing temperature. When the temperature exceeds   10 Lithosphere the liquid-point one, the heavy oil is a fluid and the velocities decrease linearly. However, as shown in Figures 14(a) and 15(a), the model overestimates the S-wave velocity of oil sands of samples #8 and #9 at the quasisolid phase, possibly because when the medium becomes solid at low temperatures, our theory does properly model the shear modulus. Figures 14(b)-16(b) show the velocities as a function of the normalized temperature. When T n is less than the liquid-point one (0.83), there is a linear trend. On the other hand, when T n is higher than the liquid-point one, the velocities deviate from this linear trend. In addition, the results also show that the porosity have no effect on the threshold of the liquid point, because the temperature has a dominant influence on the properties of the pore fluid and the effect of the skeleton is weak. However, by comparing Figure 11 (glycerol) and Figure 14(b) (heavy oil), we find that the viscosity of the fluid affects the threshold, which is different for different high-viscosity fluids.

Conclusions
We analyzed the temperature dependence of the P-and Swave velocities in high-viscosity fluid-saturated rocks. We propose a temperature-dependent CPA model for doubleporosity media, which incorporates Maxwell viscoelasticity and the David-Zimmerman models. This model effectively describes the temperature-dependent elastic properties of rocks saturated with heavy oil. In particular, the model takes into account the effects of cracks on the velocities varying with temperature.
The results show that with increasing temperature, the inelastic dispersion of the rock moduli decreases sharply,   11 Lithosphere when the high-viscosity fluid approaches a fluid phase, and the contribution of its shear modulus is negligible. The Pand S-wave velocities of four rocks saturated with a highviscosity fluid, as a function of temperature, are properly predicted with the model. The results show that the wave velocities decrease with increased temperature. It is found that the velocities varying with the temperature can be divided into two regimes according to the fluid viscosity threshold (liquid point). One regime has a linear trend and the other a nonlinear one. The porosity has no effect on the threshold, and the behavior is mainly dictated by the viscosity of the fluid.
Here, we consider the frequency-dependent characteristics of the ultrasonic wave velocities in oil sands and glycerol-saturated tight sandstones, but the method can be extended to other lithologies, e.g., carbonate rocks.

Data Availability
The data of the modeling results can be accessed by contacting the corresponding author.

Conflicts of Interest
The authors declare that they have no conflicts of interest.