Sensitivity of seismic properties to temperature variations in a geothermal reservoir

Geophysical characterization plays a key role for the definition of the deep structures of geothermal reservoirs and the consequent assessment and validation of the geothermal conceptual model. Seismic methods may provide a valuable contribution for this purpose. This involves a deep and reliable understanding of the sensitivity of seismic-wave propagation to physical and temperature variations, with complex interactions. We present the theory and sensitivity analysis based on rock's mechanical Burgers model including Arrhenius temperature equations, integrated with Gassmann model for fluid saturated porous rocks, pressure effects for bulk and shear moduli, as well as permeability and squirt flow effects. Assuming a temperature gradient model, the analysis applied at low seismic frequencies compares the interpretation of the sensitivity effects for different typical seismic elastic quantities, showing the different performance in relation to physical effects, including melting, supercritical conditions, and observability obtained in different temperature regions. With a quantification of the physical properties, the results of the study show that in deeper zones the main expected contributions in terms of variations in seismic velocity, moduli and seismic attenuation due to temperature come from melting transition, while in shallower porous fluid-saturated formations the trends are governed by pressure effects, with minor contributions of permeability and possible effects related to soft porosity. The new calculated elastic moduli are complex-valued and frequency-dependent, and temperature dependent through the fluid properties. In this complex scenario, not always the increments in the velocity and elastic wave moduli correspond to an increment in the temperature. Moreover, with mobility decreasing as a function of depth, the analysis shows that the shear quality factor is sensitive to permeability, which introduces moderate effects for velocity and attenuation of shear waves. The analysis applies to active exploration seismic and passive seismology.


Introduction
Seismic methods may provide a valuable contribution for the geophysical characterization of geothermal reservoirs, either using exploration approaches (Batini et al., 1983;Niitsuma et al., 1999) or passive seismology to image the subsurface, obtain velocity information and monitor the geothermal reservoir (e.g., Blanck et al., 2016;Majer et al., 2007). This task requires a deep and reliable understanding of the sensitivity for seismic-wave propagation to physical and temperature variations, with complex interactions of the interrelated effects. This is relevant in particular for deep-drilling projects, where supercritical fluid conditions can be encountered (Farina et al., 2016;Dobson et al., 2017;Reinsch et al., 2017) and prediction, for example by reverse VSP (RVSP) (Poletto et al., 2011;Poletto and Miranda, 2004) may play a key role.
Several works consider seismic wave propagation in hot geothermal rocks worldwide (e.g., Cermak et al., 1990;Kristinsdóttir et al., 2010;Vinciguerra et al., 2006), in the presence of temperature and fluids. Jaya et al. (2010) analysed petrophysical experiments on Icelandic geothermal rock samples at simulated in situ reservoir conditions to delineate the effect of temperature on seismic velocity and attenuation, with the goal to predict the effect of the saturating pore fluid on seismic velocity using a modified Gassman equation. In their study the temperature dependence follows solely from the thermophysical characteristic of the saturating fluid in porous rock. Iovenitti et al. (2013) and Tibuleac et al. (2013), studied the seismic-temperature distribution to test the seismic component of an exploration method calibrated by integrating geological, geophysical and geochemical experimental data, including empirical temperature -P-wave velocity relationships and sensitivity analysis after removing the effects on depth, using a geostatistical approach. More recently, seismic rheological analysis of the brittle ductile transition (BDT) and seismic propagation modeling in presence of temperature was performed by Carcione and Poletto (2013), with temperature and fluids by Carcione et al. (2014Carcione et al. ( , 2017, Farina et al. (2017), including melting and supercritical condition. The numerical algorithms developed in these studies (Carcione and Poletto, 2013;Carcione et al., 2014Carcione et al., , 2017Farina et al., 2016Farina et al., , 2017 can be used for seismic simulation in arbitrary geological media at variable geothermal conditions, including temperature and tectonic effects at depth. After understanding the seismic behavior in geothermal environments, a model-based analysis of sensitivity for the elastic quantities together with the experimental study is essential for seismic characterization. In particular, this work is part of the ongoing characterization of geothermal formations by full-waveform seismic modelling including temperature, planned and performed in the framework of the European Union Horizon 2020 GEMex Project (GEMex, 2016), for the study of high-temperature geothermal zones and geothermal systems in Mexico: for engineered geothermal system (EGS) development at Acoculco and for a super-hot resource near Los Humeros (e.g., Urban and Lermo, 2013). GEMex includes the analysis of the distribution of rock modulus of elasticity and correlation to temperature, namely: comparing the spatial distribution of rock modulus of elasticity with the temperature distribution data derived from the thermo-mechanical models with the purpose to estimate deep formation temperatures from seismic and gravity surveys (GEMex, 2016).
In this work, we present the theory and numerical sensitivity analysis based on rock's mechanical Burgers model including creep-flow by Arrhenius temperature equations, integrated with Gassmann model to account for fluid saturated porous rocks, and pressure effects for bulk modulus. The analysis includes permeability effects and squirt-flow, which may introduce unrelaxed effects at frequencies higher than the seismic frequencies. The analysis presented in the second part of the paper is applied at low seismic frequency, assuming a constant-gradient model for temperature. It compares the interpretation of the characteristic sensitivity effects for different typical seismic elastic quantities, showing the different performance in relation to physical effects, including melting and supercritical, and investigates results in different temperature regions in sample-case examples. Main results are related to interpretation of differences in sensitivity calculated with attenuation and propagation velocity, with interpretation of melting, fluid saturation and pressure effects in the sensitivity curves.
The scope of this work is to provide a first basis for the seismic sensitivity analysis with temperature by numerical simulation. The analysis is representative of the wave propagation behavior at different conditions. Carcione and Poletto (2013) observed that the Burgers model is suitable to describe the transient viscoelastic creep for arbitrary media, because there is experimental evidence that linear viscoelastic models are appropriate to describe the behavior of ductile media. Gangi (1981Gangi ( , 1983 obtained exponential functions of time using linear viscoelastic models to fit data for synthetic and natural rocksalt. Chauveau and Kaminski (2008) described the effect of transient creep on the compaction process on the basis of the Burgers' model. The viscosity can be expressed by the Arrhenius' equation, accounting for thermodynamic effects, and the constants that appear in the creep rate expressions describe the properties of a specific arbitrary material at given physical conditions. For this study, we assume isotropic materials, however anisotropy is considered in Carcione and Poletto (2013), which can be further developed for sensitivity analysis purposes. For more details on the derivation of the constitutive equations the reader may refer to previous works (Carcione and Poletto, 2013;Carcione et al., 2014Carcione et al., , 2017.

The Burgers model for brittle-ductile behavior
The constitutive equation, including both the shear viscoelastic and ductile behavior, can be described with the Burgers model as reported in Carcione and Poletto (2013) and Carcione et al. (2014). The Burgers model is a series connection of a dashpot and a Zener model ( Fig. 1) and its complex shear modulus can be written as The quantities τ σ and τ ϵ are seismic relaxation times, μ 0 is the relaxed shear modulus (see below) and η is the flow viscosity describing the ductile behavior, i = √−1 and ω = 2πf is the angular frequency. The relaxation times can be expressed as where τ 0 is a relaxation time such that ω 0 = 1/τ 0 is the center frequency of the relaxation peak and Q 0 is the minimum quality factor. The limit η → ∞ in Eq. (1) recovers the Zener kernel to describe the behavior of the brittle material, while τ σ → 0 and τ ϵ → 0 yield the Maxwell model used by Dragoni and Pondrelli (1991): μ B = μ 0 (1 − iμ 0 / ωη) −1 (e.g., Carcione, 2014). For η → 0, μ B → 0 and the medium becomes a fluid. Moreover, if ω → ∞, μ B → μ 0 τ ϵ /τ σ , where μ 0 is the relaxed (ω = 0) shear modulus of the Zener element (η =∞).
The viscosity η can be expressed by the Arrhenius equation (e.g., Carcione et al., 2006;Montesi, 2007). It is related to the steady-state creep rate ϵ˙by where σ o is the octahedral stress (e.g., Gangi, 1981Gangi, , 1983Carcione et al., 2006;Carcione and Poletto, 2013), A ∞ and n are constants, E is the activation energy, R G = 8.3144 J/mol/°K is the gas constant and T is the absolute temperature. The octahedral stress is where the σ's are the stress components in the principal system, corresponding to the vertical (v) lithostatic stress, and the maximum (H) and minimum (h) horizontal tectonic stresses. The temperature is a function of depth through the geothermal gradient G. A linear approximation is T − T 0 = z G, where z is the depth and T 0 is the temperature at the surface (z = 0).

The modified Gassmann model
Gassmann's equations are used to calculate changes in seismic velocity and elastic quantities due to different fluid saturations. In this Fig. 1. Mechanical representation of the Burgers viscoelastic model for shear deformations (e.g., Carcione, 2014). σ, ϵ, μ and η represent stress, strain, shear modulus and viscosity, respectively, where η 1 describes seismic relaxation while η is related to plastic flow and processes such as dislocation creep.
work we assume that the porous material is isotropic, and homogeneous. The Gassmann bulk and shear moduli are where ϕ is the porosity, K m and μ m are the bulk and shear moduli of the drained matrix, and K s and K f are the solid and fluid bulk moduli, respectively (e.g., Carcione, 2014).
To account for the pressure dependence, we express the dry-rock bulk moduli as where g j (p d ), j = 1, 2 defines the dependence of the moduli on the differential pressure p d = p c − p, where p c is the confining pressure, p is the pore (fluid) pressure, and K 0 and μ 0 are the bulk and shear moduli at infinite effective pressure and η =∞ (or ω =∞). Using μ B in (8) means that the Burgers shear viscosity is included. The simplest form of function g, in good agreement with experimental data, is (Kaselow and Shapiro, 2004), where a j and p* j are parameters. It is g j = 1 for p d → ∞ (e.g., very high confining pressure) and g j = a j for p d → 0 (pore pressure equal to the confining pressure). The bulk density is where ρ s and ρ f are the grain and fluid densities, respectively. In the following analysis we distinguish between stiff grain porosity and soft compliant porosity in the gap area of grain contact. The compliant porosity is typically so smallnearly 0.001 for most rocksthat the total porosity ϕ can be assumed to be equal to the stiff porosity.

Phase velocity, attenuation and wave modulus
The phase velocity and attenuation, or dissipation factor (inverse of the quality factor), including the Burgers, Biot, permeability and squirtflow losses (see next sections), are and where v c is the complex velocity (e.g., Carcione, 2014). For shear waves and , where is the rock tortuosity and κ is the permeability. In our simulations we assume = − − ϕ 1 0.5(1 1/ ) (Mavko et al., 2009;Berryman, 1980). Eq. (14) can be reformulated as where the quantity = p T κ η ( , ) ( / ) f is mobility, ratio of permeability and viscosity (Batzle et al., 2006), introducing dispersion and attenuation effects in the shear-wave Eqs. (12) and (13), as well as in Eq. (18) for compressional waves, depending on pressure and temperature through fluid properties and permeability. Manning and Ingebritsen (1999) inferred permeability from thermal modeling and metamorphic systems suggesting the following dependence with depth z, where T 0 = T(0) is the surface temperature, z is the depth in km and the permeability is given in m 2 . With constant-gradient approximation, the second expression assumes a linear geothermal law, T − T 0 = z G. The complex velocity of the P waves is obtained from the following second-order equation in v c 2 : (2 ) , (e.g., Carcione, 2014 Eq. (7.324)). The P-wave and S-wave stiffness moduli are given by respectively.

Sensitivity analysis
The approach we adopt to investigate the sensitivity of the seismic properties to temperature T is as follows. Consider the quantities where ρ, are formation density, and compressional and shear velocities, respectively, and Q is the quality factor (12) accounting for attenuation related to temperature. We define the stiffness (modulus, denoted by subscript 'M') and impedance (denoted by subscript 'I') sensitivities as and that of attenuation by where the subscript and superscript (J = P, S) denote the compressional-or shear-wave type index. Eq. (21) expresses the effect of temperature on the stress-strain relations through the moduli of elasticity, while Eq. (22) refers to the radiation impedance, since it includes density and not only wave velocity. Note that the sensitivity for the Bulk modulus is given by The quantities in Eqs. (21) and (22) can simply be expressed as a function of ρ, (∂ρ/∂T), v J and ∂ ∂ v T ( / ) J . In poro-viscoelastic media, for a given type of rock and saturating fluid, it is in general where ϕ is porosity, ρ f is the fluid density and p is pressure. In general at variable depths and in the proximity of melting conditions we also consider the dependence ϕ = ϕ(T, p). Moreover, the quantities v J and Q J also depend on local stress conditions and typically they exhibit non negligible dispersion effects, since they depend on frequency ω. Namely, we have = = T ω σ ϕ p κ Q Q T ω σ ϕ p κ ν ν ( , , , , , ), ( , , , , , ), where we introduced the octahedral stress σ o and the permeability κ is included. With these premises, we analyze the sensitivity quantities s K J ( ) where the subscript K = M, I, Q. These curves provide tools to evaluate the reliability and the effectiveness of the temperature analysis methods by seismic signals. To achieve this task using different seismic quantities, we calculate and compare also the relative-sensitivity curves T ( ) for the investigated quantities, generically denoted by W J , with respect to the value of the same quantity at the same temperature, which provides an estimator of the performance of the different quantities W J to characterize seismically the geothermal variation. For example, with slow and negligible variations of the impedance = W T ρv ( ) J J in limited temperature intervals, Eq. (27) gives an approximation of twice the magnitude of the acoustic reflection coefficient per unit temperature increment, and can be used to evaluate the impact of the investigation by the seismic reflection response once the dependence of the seismic quantity on temperature and the local temperature gradient are known, hence for inversion purposes.
Finally, note that the reciprocal of the sensitivity of seismic quantities to temperature, say, of a given measurable quantity W J , provides the rate of variation (sensitivity) of temperature with respect to said measurable seismic quantity, which can be investigated and utilized for evaluation of stability conditions in the seismic prediction of geothermal temperature distribution and variation between depth intervals. When the reciprocal sensitivity of Eq. (28) is low, a variation in the seismic quantity correspond to a lower variation in temperature, and the prediction is locally more stable.

Physics of the fluid-saturated rocks
The rheological conditions we study to analyze seismically geothermal fields include solid rock properties, fluid properties, temperature, pressure, tectonic conditions, porosity and permeability versus depth and temperature. Possible squirt flow effects are also evaluated.

Pressure, tectonic stress and thermal parameters
The pressure dependence at seismic frequencies (Eq. where p d is given in MPa  and the constants a j and p* j in Eq. (9) are calculated from Popp and Kern (1994). The shear seismic loss parameter is obtained from empirical equations derived by Castro et al. (2008) for the crust in Southern Italy. They report Q 0 = 18.8 f 1.7 for the upper crust and up to a frequency of 10 Hz. In the examples we consider a frequency of f = 3 Hz, with ω 0 = 2 πf, which gives Q 0 = 122. The temperature is a function of depth through the geothermal gradient G as T = (20 + z G × 10 −3 ), where z (m) is depth and G = 50, 60 and 90°C/km in our calculations for different examples. The lithostatic stress is = − = − σ ρ gz p v c , where ρ = 2400 kg/m 3 is the average density and g = 9.81 m/s 2 is the gravity constant. To obtain the octahedral stress (4) we consider a simple model based on the gravity contribution at depth z. The horizontal stresses are estimated as where ν = ν(K 0 , μ 0 ) is the Poisson ratio at infinite effective pressure. The factor ν/(1 − ν) lies between 0.25 and 1 for ν ranging from 0.2 to 0.5, with the latter value corresponding to a liquid (hydrostatic stress). The parameter ξ ≤ 1 has been introduced to model additional effects due to tectonic activity (anisotropic tectonic stress) (Carcione and Poletto, 2013). Furthermore, we consider A ∞ = 100 (MPa) −n s −1 , E = 134 kJ/mol and n = 2, and take ξ = 0.8. The above degree of stress anisotropy is consistent with values at prospective depths provided by Hegret (1987) for the Canadian Shield, and in agreement with data reported in Engelder (1993, p. 91).

Fluid physical properties
Without loss of generality, in our examples the geothermal fluid is pure water. The water properties as a function of pressure and temperature are obtained from the fluid thermo-physical database provided in the website of the National Institute of Standards and Technology (NIST), collected from laboratory measurements by Lemmon et al. (2005). In "Thermophysical Properties of Fluid Systems", we choose water (1)  With only liquid, a state of hydrostatic pore pressure is given by = p ρ gz f , where ρ f = 1000 kg/m 3 is an average fluid density. In the presence of different, liquid, vapour and supercritical phases at depth, we use an iterative method to calculate the hydrostatic pressure by NIST (Farina et al., 2016). In Fig. 2 we consider a depth range [5,15] km, where pore pressure and temperature vary from 50 to 150 MPa and 300 to 900°C, respectively (in this example the geothermal gradient is 60°C/km). The experimental density, sound velocity and viscosity of water are shown together with the pressure and temperature profiles. Compare these values to the ones at ambient conditions, defined by a temperature of 20°C and a pressure of 0.1 MPa: a water density of 998 kg/m 3 and a sound velocity of 1482 m/s.

Porosity and permeability
Let us consider now variations in the rock porosity ϕ, assuming stiff porosity approximation ϕ ≅ ϕ s . This implies that the bulk and shear moduli depend on porosity by Gassmann model as well as on the permeability. The dry bulk and shear moduli of the samples are determined by the Krief model (Krief et al., 1990), where K s and μ s are the bulk and shear moduli of the solid. Permeability is obtained as Poletto et al. Geothermics 76 (2018) 149-163 (Carcione et al., 2000), where R s is the average radius of the grains.
Here, we assume R s = 20 μm. As we can see below, this relation is in agreement with rheological estimations with temperature in different types of rocks, and we assume it as a generalized approximated relation, with possible deviations for metamorphic rocks. Moreover, significant deviations can be expected for permeability in formations with secondary crack porosity but low stiff porosity in rocks with flow paths in faulting and fractures, both when considering natural systems (Hickman et al., 1995;Ito and Zoback, 2000) and EGS systems with enhanced fracturation (Majer et al., 2007;Hashida et al., 2001). To account for the dependence of permeability from depth and temperature, we invert Eq. (32) for porosity by solving where = κ κ 45 and κ = κ(T) is given from Eq. (16). The plot of permeability as a function of porosity shown in Fig. 3a is in agreements with the curves estimated empirically for sandstone and carbonate rocks by Ehrenberg and Nadeau (2005, in particular in Fig. 4 of their paper). In this example, the dependence on temperature is calculated as T = zG + T 0 with T 0 = 20°C and a gradient G = 50°C/km in order to avoid melting zones. In order to display the curves in the same numerical range, Fig. 3b shows the permeability expressed in (mDarcy × 10) unit and porosity (%) versus temperature in the depth interval [2,8] km. In this analysis, according with Ehrenberg and Nadeau (2005), we neglect the shallower layers with higher porosity and, in agreement with Fig. 3b, we consider examples with the porosity in the range ϕ = [0, 5] (%).

Stiff and soft pores and squirt-flow effects
The squirt flow interaction model between stiff and soft pores takes into account the fact that the pore space of many rocks has a binary structure composed of relatively stiff pores, which constitute the majority of the pore space, and relatively compliant (or soft) pores, which are responsible for the pressure dependency of the poroelastic moduli. Fluid saturates both stiff and soft pores. When the frequency is higher than the so-called characteristic squirt relaxation frequency f CS , the fluid pressure does not have enough time to equilibrate between stiff and compliant pores during a half-wave cycle. Above f CS the system is in the so called unrelaxed state. Then, compliant pores at the grain  contacts are effectively isolated from the stiff pores and hence become stiffer with respect to normal (but not tangential) deformation.
In order to model the frequency dependency of the partially-relaxed moduli, Gurevich et al. (2010) assumed a geometrical configuration by which a compliant pore forms a disk-shaped gap between two grains, and its edge opens into a toroidal stiff pore (Fig. 4). Gurevich et al. (2009Gurevich et al. ( , 2010 analyzed the ultrasonic behaviour, and the low-, intermediate-and high-frequency approximations for squirt-induced attenuation. They obtained the modified partially-relaxed (at sufficient low frequency) dry moduli, i.e., whereby soft pores are fluid-filled whereas stiff pores are dry, as where K m and μ m represent the dry-rock bulk and shear moduli at the confining pressure p c , K h is the dry-rock bulk modulus at a confining pressure where all the compliant pores are closed, i.e., that of a hypothetical rock without the soft porosity, and ϕ c is the compliant (soft) porosity. For a more detailed description of the numerical modeling approach see Carcione and Gurevich (2011). The key quantity in Eqs.
and k is the wavenumber, R is the radius of the crack and h is its thickness (Fig. 4). When the fluid modulus satisfies is an effective viscosity. The peak relaxation frequency of the squirtflow model is using the approximations K h ≈ K m and assuming ≫ K ϕ ( ) s c (Carcione and Gurevich, 2011). Hence, the peak frequency decreases with increasing fluid viscosity and decreasing aspect ratio (h/R) of the crack.
Using Eqs. (34) in the Gassmann model Eqs. (5)- (7) gives the modified squirt-Gassmann moduli. The explicit functional form of α and M on K m is in fact convenient for replacing K m by the modified matrix (or frame) complex modulus K including the squirt-flow mechanism (Eq. (34)). Similarly, μ m is replaced by μ. The new moduli are complexvalued and frequency-dependent, and, relevant for our study, also temperature dependent through the fluid properties.
As discussed by Carcione et al. (2018a,b), the squirt-flow model is consistent with Gassmann's theory in the low-frequency limit, and with Mavko-Jizba unrelaxed moduli in the high-frequency limit (Mavko and Jizba, 1991). All the parameters of the model have a clear physical meaning. There is only one adjustable parameter: the aspect ratio of compliant pores (grain contacts) h/R. However, the model approximations in different frequency regions are different for different fluid phases, i.e., not only fluid but also gas (Gurevich et al., 2010;Carcione and Gurevich, 2011) and supercritical. The squirt physical effect, here introduced for a preliminary evaluation for the purposes of the sensitivity analysis, needs further investigations with multi-phase fluids to evaluate its relevance at seismic frequencies with temperature.

Case study for seismic and physical quantities
Using rock and geothermal parameters of the reference literature, we calculate the seismic elastic quantities to obtain characteristic sensitivity curves for a small set of explanatory physical models. We compare the sensitivities (∂ρ/∂T), ∂ ∂ v T ( / ) J and (∂Q J /∂T) together with the normalized sensitivity curves of Eqs. (21) and (22), for a uniform formation with background temperature-unperturbed compressional velocity = v 6670 P m/s, shear velocity = v 3851 S m/s and unperturbed density ρ = 3000 kg/m 3 , for the following parametrizations in the lowfrequency approximation at seismic frequencies (Table 1). Burgers and thermal parameters are listed in Table 2. The discussion is focused on sensitivity calculated by Burgers model including shear loss with temperature, Gassmann model with fluid saturation, as well as bulk modulus dependence on pressure (as expressed by Eq. (9)), stiff porosity and permeability. Examples with squirt flow are presented. Fig. 5a shows the fluid density calculated with the values of Tables 1 and 2 and the temperature gradient of 90°C/km shown in Fig. 6a. Fig. 5b plots the water pressure versus temperature, where the supercritical zone for temperature and pressure in excess to T = 374°C and 22.1 MPa, respectively, is evidenced. Note that in the region of (a) corresponding to the supercritical zone in (b), the density is lower and consequently the pressure increases with a slower trend as a function of depth and temperature in the supercritical zone of (b).

Seismic velocity
Signals are calculated at the reference frequency f = 10 Hz (Table 2), using the linear temperature-depth model T(z) = z G × 10 −3 + T 0 , where T 0 = T(0) = 20°C, with constant gradient G = 90°C/km of Fig. 6a. Fig. 6b shows the wave velocity for P and S waves. Without melting (i.e., neglecting the Burgers viscosity), the wave velocities have minor variations.
In the following figures we interpret contributions due to different physical effects in the plots of the seismic quantities versus temperature. First we compare different behaviors in the responses obtained by velocity and elastic moduli. Fig. 7 represents v P and v S in the presence and absence of saturating geothermal fluids, therefore with and without porosity and with and without pressure effects for the moduli. Here and in the following we use the term 'without porosity' or 'zero porosity' to denote negligible porosity, e.g., less than 0.1 %. In Fig. 7a and b we Fig. 4. Sketch of the squirt-flow model, where two sandstone grains in contact are shown. The soft pores are the grain contacts and the stiff pores constitute the main porosity. The quantity R is the radius of the disk-shaped soft pore (half disk is represented in the plot) (modified after Gurevich et al., 2010). observe the large effect for v P and v S , respectively, in the melting zone, similar for all the curves of the same panels. In the curve of v P calculated with porosity, shown by the red line in Fig. 7a (stiff porosity ϕ = 0.05), we observe the zone where the fluid saturation effect, indicated as the 'Gassman zone', is more evident, while the curve without porosity is flat in that zone. From this result we deduce that, depending on porosity value, we can use velocity information to investigate the geothermal fluid saturation effects in the model, also when the pressure effects for the bulk modulus are not included. The curve with pressure effects for the bulk modulus (dashed line) presents variations also at lower temperatures, below the melting zones. This corresponds to different sensitivity curves, as we will see what follows. Fig. 7b shows the corresponding curves for the S-wave components, where the Gassmann effect is much less evident. In this case, different from P-wave, the shear velocity tends to zero at high temperatures as expected beyond the melting zone, where the rock is fluid.

Seismic stiffness and density
The corresponding curves calculated with the P-wave and S-wave elastic moduli show similar trends with melting. However, as it can be observed in Figs. 8a and b, the same considerations in relation to Gassmann fluid saturation effects in the presence of porosity are not valid for the sensitivity analysis with the elastic moduli. In Fig. 8a we show the P-wave elastic modulus E P with and without porosity. and pressure effects. We observe that the presence of fluid changes the Pwave elastic modulus in the curves without pressure effects, which becomes lower, but these curves are parallel, therefore they present the same sensitivity. Moreover in the absence of the pressure correction they are both flat below the melting zone. Fig. 8b shows the S-wave elastic modulus E S with and without porosity and pressure effects for the bulk and shear moduli. In this case the curves with and without porosity are superimposed.
In the presence of porosity and fluid saturation the velocity changes and also the density changes, in such a way that the change in v P S , 2 is inversely proportional to that of density, and this creates a compensation effect in the elastic moduli. The compensation effect observed in the flat regions of Fig. 8 is explained using Fig. 9. The result is that there is not variation in E P,S relative to temperature for Gassmann effects. Fig. 9 shows these trends in normalized curves calculated with porosity and without pressure, in this case using both P-wave and S-wave elastic moduli. We compare v P 2 , v S 2 and 1/ρ by amplitudes normalized at the temperature origin (with unit relative amplitude at T = 0). In the region below the melting zone these curves are superimposed, since the product of density and v 2 eliminates the variations relative to temperature due to fluid saturation.

Temperature as a function of seismic quantities
In temperature-velocity regression analysis, e.g., such as in Iovenitti et al. (2013), it is sometimes convenient to exchange the plot axes, and to represent the temperature versus velocity, or other seismic parameters. Fig. 10a shows the plots of temperature as a function of velocity, both using P-wave and S-wave velocities calculated with fluid saturation (ϕ = 5%) and pressure effects for the Bulk and shear wave's moduli. Fig. 10b shows the similar plots for the P-wave and S-wave moduli. Obviously, this type of representation depends on the temperature trend versus depth, in this case a gradient. Fig. 11 shows a detail of the same curves of Fig. 10 in a restricted temperature region [100, 600]°C, to better evidence the trends below the melting zone.

Permeability, mobility and attenuation
We extend the poro-viscoelastic model to include the permeability (Carcione et al., 2018a,b). Permeability, according with the depth dependence given by Eq. (16) decreases versus depth, and it can be assumed negligible in the melting zone for the purposes of our sensitivity analysis. The effect is governed by fluid mobility and is dispersive. Fig. 12 shows (a) the fluid viscosity versus temperature and (b) the mobility in the temperature interval [100, 800]°C, below the lower limit of the melting zone. At higher temperatures, mobility is close to zero because the permeability decreases with depth and vanishes in the proximity of the melting zone (Fig. 3b).
The analysis shows that the shear quality factor is sensitive to permeability. Fig. 13a shows the shear-wave Q S quality factor calculated in the temperature interval [0, 800]°C without and with permeability. The signal frequency is 200 Hz. In this example the unperturbed intrinsic attenuation is low. For permeability, we use the variable function of depth given by Eq. (16). The result shows observable variations with respect to the case without permeability, especially at lower temperatures, where permeability is higher. To evidence possible effects also at higher depths, we test also the approach using constant permeability κ C = 1.5 × 10 −14 m 2 (blue curve in Fig. 13a). In this case Fig. 7. Plot of (a) P-wave velocity v P , calculated with and without porosity, and pressure effects for the bulk and shear moduli. In the curve with porosity we observe the fluid Gassmann effect, while the curve without porosity is flat in the Gassmann zone. (b) S-wave velocity v S , calculated with and without porosity, and pressure effects for the bulk modulus. The curves with and without porosity and both without pressure are superimposed. (See Table 3 for the list of the symbols).  8. (a) P-wave modulus E P curves, calculated with and without porosity, and pressure effects for the bulk and shear moduli. The curves without pressure are parallel, hence they present the same sensitivity to temperature variation. (b) S-wave modulus E S curves, calculated with and without porosity, and pressure effects for the bulk and shear moduli.

Squirt flow
Assuming the presence of stiff and soft porosity (Fig. 4), the squirt flow modulus is calculated using the following parameters: solid density ρ s = 3000 kg/m 3 and K s = 74.2 GPa (this value is deduced from Table II of Popp and Kern (1994), intrinsic velocity data at 200 MPa), h/ R = 0.00001, K h = 66.2 GPa and ϕ c = 0.00001. Closure of cracks with confining pressure is reflected in the values of the compliant porosity given in Table II of Popp and Kern (1994), ranging from 0.28% at 12 MPa to 0.01% at 200 MPa. We investigate possible dispersion effects introduced by squirt flow. We observe that K h , a key value that is the bulk modulus of the hypothetical rock without compliant porosity (Gurevich et al., 2010), determines also the trends of K at low ϕ c in Eq. (34). To prevent from distorted physical effects at low pressure, otherwise we obtain large differences in K and K m at shallower depths and lower temperatures with negligible ϕ c , we have to consider the variability in the compliant porosity with depth. Following Gurevich et al. (2010) where for K f we try both K f and K * f . We also take advantage from the laboratory results reported by Popp and Kern to infer a decay curve for the compliant porosity. These curves are compared in Fig. 14 in the temperature range [0, 900] (°C). Fig. 15 shows the velocity of (a) P-waves and (b) S-waves calculated with only Burgers (B) at fixed frequency 10 Hz, and Burgers plus squirt flow calculated with variable ϕ c estimated by Pop and Kern (1994) (P & K), and by Eq. (38) using K * f . The model includes porosity and pressure effects. More evident for the P-waves (a), below the melting zone (in this case this is the interpretation zone) the velocity calculated with only Burgers increases when the Burgers model is used together with squirt flow. Fig. 16 shows the corresponding quality factor Q P and Q S of (a) P-waves and (b) S-waves, respectively, calculated with only Burgers (B) at fixed frequency 10 Hz, and Burgers plus squirt flow calculated with variable ϕ c estimated by Pop and Kern (1994) (P & K), and by Eq. (38) using K * f . Although the magnitude of the simulated velocities and attenuation can be revised and could be matter of further evaluation for a suitable choice of the rheological parameters in order to calibrate the model in a geothermal context, this result shows that the squirt flow may introduce effects at low frequencies (in this case 10 Hz) depending on compliant porosity estimate.   Fig. 10b, showing the curves in the temperature interval [100, 600]°C, below melting. The curves are calculated for both for P-wave (solid line) and S-wave (dashed line) moduli E P and E S .

Sensitivity analysis and interpretation
The numerical sensitivity analysis is performed using the Burgers-Gassmann poro-viscoelastic model with pressure, for four different conditions summarized in Table 3: with and without porosity effects (P05 and P00, respectively), and without and with pressure effects (G1 and G, respectively) for the bulk and shear moduli. Namely, G1 and G denote that no pressure effect is accounted for, and that the pressure effect is accounted for, respectively. For example, the notation P05G in the figures denotes that the porosity ϕ = 0.05 is used to account for the Gassmann behavior and pressure effects are accounted for to calculate the bulk modulus (see Table 3). Each panel compares these four different sensitivity curves for a selected quantity, by showing the absolute value of the sensitivity in Figs. from 17 to 21.
The interpretative analysis shows the different temperature regions where the sensitivity variations related to the physical effects are more important, for the different model curves. Not all the curves in the panels are affected by the same effects. In all the panels, the Burgers melting is the more important effect. Superimposed are the Gassmann fluid effect with porosity and fluid saturation, and the effect related to pressure-induced bulk modulus variation. Permeability introduces moderate effects for velocity and attenuation of shear waves. Fig. 17a shows the sensitivity curves for v P . In the presence of porosity and fluid-saturation we observe Gassmann effects. With pressure-bulk modulus correction we observe evidence of trends at low temperatures. The arrows indicate schematically the zones in which the different effects are more relevant for the sensitivity, displayed by absolute value. Similar to Fig. 17a, Fig. 17b shows the sensitivity curves for v S . Note the different extension along the temperature axis of the Burgers melting sensitivity region with respect to v P . Fig. 18a shows the sensitivity curves for the P-wave elastic modulus E P . A weaker Gassmann effect is observable, only with both porosity and pressurebulk correction. Fig. 18b shows the sensitivity curves for the S-wave elastic modulus E S . As expected, Gassmann effects are not observable in the shear sensitivity plots.
The analysis is also applied to attenuation effects. Fig. 19a shows the sensitivity curves for the Q P factor. Prevalent effect of Burgers melting can be observed, also at high temperatures, where the sensitivity of Q P increases. Fig. 19b shows the sensitivity curves for Q S . All the curves are superimposed, and only the effect of Burgers melting is present, in this case only in the melting zone around the peak approximately at 500°C. For Q S , we observe that there is no increase in the sensitivity for increasing temperature as for Q p , since after melting the shear waves do not propagate in the magma fluid. Finally, in Fig. 20a we see the sensitivity (absolute value) of the shear quality factor Q S calculated without and withvariable and constantpermeability. These sensitivity data are calculated using the signals shown in Fig. 13a in the Fig. 12. Plot of (a) fluid viscosity and (b) mobility in the temperature interval [100, 800]°C, below and at the lower limit of the melting zone. At higher temperatures mobility is close to zero because permeability decreases with depth and vanishes in the proximity of the melting zone.  temperature interval [0, 800]°C. We observe, especially at shallower depths and lower temperatures, significant variations of sensitivity for the curves calculated with variable permeability, decreasing with depth, and constant permeability (see the example of Fig. 13a). Fig. 20b shows the sensitivity curves calculated for density ρ. Relevant variation of its sensitivity calculated with porosity and fluid saturation is observable, and interpreted also as related to supercritical effects assuming as geothermal fluid pure water.
All these examples show that the different quantities provide, with different extents and case by case, better estimations of sensitivity in different temperature regions. In Fig. 21 we compare the absolute values of the normalized sensitivity curves of different elastic quantities. The curves in this superposition cover with different responses different temperature regions. Note that, as observable in the previous plots, the peaks of the sensitivity are at different temperatures for velocity (at approximately 800°C) and quality factor (at approximately 500°C). We may better observe the out-of phase trends of attenuation and dispersion in the next figure.
In Fig. 22a we compare the corresponding relative sensitivity responses ( by Eq. (27)) for the same quantities of Fig. 21, i.e., each curve represents its relative variation with respect to its physical value at the given temperature per temperature degree. In this case the sensitivity curves are plotted with positive and negative signs, to show the polarity of the relative variations. In this figure, the out-of-phase behavior of velocity and Q-factor is more evident. This is similar to the fact that for causal physical signals dispersion and attenuation are Kramers-Krönig pairs (e.g., Sun et al., 2009). We may see that a relevant   observability effect is obtained in the melting zone using the Q-factor, and in general more prominent with shear components. Fig. 22b shows the estimated, and approximated as previously discussed, reflection coefficient calculated by × T Δ J at each T value using Eq. (27) for compressional-and shear-wave impedances and a temperature interval ΔT = 10°C, as with a step ΔT between two uniform temperature zones. Also in this plot we display the curves with positive and negative signs. This provides an estimation of the reflection response related only to the temperature model.
Finally, it is typically convenient using the reciprocal of sensitivity to predict (in stable regions, i.e., where the sensitivity is different from zero) temperature variations for an increment of velocity. Fig. 23 shows the predicted temperature variation     S m/s. In both the P-wave and S-wave plots we observe zones where an increment in the velocity can correspond both to an increment or to a decrease of the temperature in the different regions interpreted in the figures. Note that higher sensitivity means higher detectability of velocity changes induced by temperature, and a more stable result for a given velocity variation when we predict temperature from velocity. For example, a large P-wave velocity variation of 100 m/ s at = v 5800 P m/s in the melting zone in Fig. 23a corresponds to a decrease in the temperature of approximately 10°C, while in the pressure zone it corresponds to an increase in temperature of approximately 50°C.

Discussion and research perspectives
The paper presents a comprehensive analysis of temperature effects for fluid-saturated rock together with a theoretical basis and models for the seismic characterization of geothermal formations. The interpretation of results points out different trends and effects in the sensitivity analysis. These are related to different models corresponding to specific physical effects. The interaction of these physical conditions and effects is typically complex. In this analysis, the choice and definition of the temperature distribution map, approximated by a constant-gradient model for our purposes, is of great importance.
The characteristic sensitivity examples shown here for a case study are numerically calculated at fixed parameters, using a set of physical configurations, and low frequency, and are not exhaustive for a characterization of geothermal systems belonging to different and much more complex geological scenarios. For example, the change of the rock type and of its Arrhenius parameters, as well as tectonic stresses, may change the melting temperature and this may cause a different superposition of the physical effects in the sensitivity curves versus temperature. The change of the geothermal fluid properties (e.g., Jaya et al., 2010) may change the supercritical point, here assumed to be that of pure water, hence pressure and density curves used in the numerical calculations. Further corrections and improvements can be introduced to take into account further physical effects and relationships, as well as experimental evidences in the modeling and sensitivity analysis of geothermal seismic properties and wave fields, such as those of the Mexican high-enthalpy areas (GEMex Project, 2016). In this case, the main targets will be to characterize seismically the super-hot geothermal systems, investigate the permeability and fracturation conditions, evaluate the possible presence of fluids at supercritical conditions, and contribute to map possible magmatic zones (BDT) interpreted in the proximity of the investigated areas, thus supporting geothermal exploration and future exploration and production drilling.

Conclusions
Understanding the sensitivity of seismic quantities to temperature is of great importance for the seismic characterization of geothermal reservoirs. Especially at high temperatures, detection and monitoring of melting and supercritical zones, as well as influence of pressure on the bulk and shear moduli require appropriate sensitivity analysis. In this paper we present the Burger-Gassmann theory following previous studies and numerical-code developments, including permeability and involving squirt-flow effects to some extent, and study characteristic sensitivity curves in the low-frequency approximation. Results show the  = v Δ 58 S m/s in the S-velocity interval [500, 3500] m/s. In both the compressional and shear panels we observe zones where an increment in the velocity can correspond both to an increment or to a decrease of the temperature. different observability by different elastic components, with different prevalence of the physical effects in different temperature regions. This suggests the use of an integrated analysis by more seismic elastic quantities for the characterization of geothermal areas, which can be applied either to exploration or to passive seismology data, including volcanic environments.
The characteristic sensitivity is here calculated for a set of physical models. Based on a quantification of the physical properties, the results show that in deeper zones the main expected contributions in terms of variations in seismic velocity, moduli and seismic attenuation due to temperature come from melting transition, while in shallower porous fluid-saturated formations the trends are mainly governed by pressure effects, with minor contributions of permeability and possible effects related to the compliant soft porosity. In the region corresponding to the supercritical zone, the fluid density is lower and consequently the pressure increases with a slower trend as a function of depth and temperature. Without melting (i.e., neglecting the Burgers viscosity), the wave velocities have minor variations. Depending on porosity, we can use velocity information to retrieve the fluid saturation. The trend including pressure effects in the bulk and shear moduli presents variations even at low temperatures. The Gassmann effect is less evident in the S-wave velocity, which tends to zero at high temperatures due to melting, as expected. In the curves calculated without pressure effects for the bulk and shear moduli, the presence of fluid changes the P-wave elastic modulus which becomes lower than that calculated in the absence of fluid, but these curves are parallel, therefore they present the same sensitivity. For the S-wave elastic modulus with and without porosity and pressure effects the curves are practically superimposed. In the presence of porosity with fluid saturation the velocity and the density change, in such a way that the P-wave modulus is almost constant with temperature regarding the Gassmann effects. In the analysis of temperature as a function of seismic quantities by reciprocal sensitivity, not always the increments in the velocity and elastic wave moduli correspond to an increment in the temperature. For example, the same increment in the S-wave velocity may correspond to an increase in the temperature in a zone where pressure effects are observed and to a decrease in the temperature in the melting zone.
The fluid viscosity decreases initially as a function of temperature and then increases slowly in the supercritical zone. At high temperatures, the fluid mobility is close to zero because the permeability decreases with depth and vanishes in the melting zone. The analysis shows that the shear quality factor is sensitive to permeability. Permeability introduces moderate effects for velocity and attenuation of shear waves. We observe these effects, especially at shallower depths and low temperatures, for the curves calculated with variable permeability, decreasing with depth. Moreover, assuming a constant-permeability model, we study the potential permeability effects for deeper zones.
In this analysis, the choice and definition of the temperature distribution map, approximated by a constant gradient for our purposes, is important. The change of the rock type and its Arrhenius parameters, as well as the tectonic stresses, may change the melting temperature and this may cause a different distribution of the physical effects, partially superimposed in the sensitivity curves. The change of the geothermal fluid properties affects the supercritical point, here assumed to be that of pure water, hence the pressure and density curves used in the calculations.
Next, we plan to apply the analysis to real cases, such those of the Mexican high enthalpy regions, where the main targets are to characterize seismically the super-hot geothermal systems, including the temperature, evaluate the possible presence of supercritical-fluid conditions, and contribute to map possible magmatic zones interpreted in the proximity of the investigated areas.