Numerical modeling of fluid effects on seismic properties of fractured magmatic geothermal reservoirs

Seismic investigations of geothermal reservoirs over the last 20 years have sought to interpret the resulting tomograms and reflection images in terms of the degree of reservoir fracturing and fluid content. Since the former provides the pathways and the latter acts as the medium for transporting geothermal energy, such information is needed to evaluate the quality of the reservoir. In conventional rock physics-based interpretations, this hydro-mechanical information is approximated from seismic velocities computed at the low-frequency (field-based) and high-frequency (labbased) limits. In this paper, we demonstrate how seismic properties of fluid-filled, fractured reservoirs can be modeled over the full frequency spectrum using a numerical simulation technique which has become popular in recent years. This technique is based on Biot’s theory of poroelasticity and enables the modeling of the seismic velocity dispersion and the frequency dependent seismic attenuation due to waveinduced fluid flow. These properties are sensitive to key parameters such as the hydraulic permeability of fractures as well as the compressibility and viscosity of the pore fluids. Applying the poroelastic modeling technique to the specific case of a magmatic geothermal system under stress due to the weight of the overlying rocks requires careful parameterization of the model. This includes consideration of the diversity of rock types occurring in the magmatic system and examination of the confining-pressure dependency of each input parameter. After the evaluation of all input parameters, we use our modeling technique to determine the seismic attenuation factors and phase velocities of a rock containing a complex interconnected fracture network, whose geometry is based on a fractured geothermal reservoir in Iceland. Our results indicate that in a magmatic geothermal reservoir the overall seismic velocity structure mainly reflects the lithological heterogeneity of the system, whereas indicators for reservoir permeability and fluid content are deducible from the magnitude of seismic attenuation and the critical frequency at which the peak of attenuation and maximum velocity dispersion occur. The study demonstrates how numerical modeling provides a valuable tool to overcome interpretation ambiguity and to gain a better understanding of the hydrology of geothermal systems, which are embedded in a highly heterogeneous host medium.

pathways and the latter acts as the medium for transporting geothermal energy, such information is needed to evaluate the quality of the reservoir. In conventional rock physics-based interpretations, this hydro-mechanical information is approximated from seismic velocities computed at the low frequency (field-based) and high frequency (lab-based) limits. In this paper, we 5 demonstrate how seismic properties of fluid-filled, fractured reservoirs can be modeled over the full frequency spectrum using a numerical simulation technique which has become popular in recent years. It is based on Biot's theory of poroelasticity and enables the modeling of the seismic velocity dispersion and the frequency dependent seismic attenuation due to waveinduced fluid flow. These properties are sensitive to key parameters such as the hydraulic permeability of fractures as well as the compressibility and viscosity of the pore fluids. Applying the poroelastic modeling technique to the specific case of a 10 magmatic geothermal system under stress due to the weight of the overlying rocks requires careful parameterization of the model. This includes consideration of the diversity of rock types occurring in the magmatic system and examination of the confining pressure-dependency of each input parameter. After the evaluation of all input parameters, we use our modeling technique to determine the seismic attenuation factors and phase velocities of a rock containing a complex interconnected fracture network, whose geometry is based on a fractured geothermal reservoir in Iceland. Our results indicate that in a magmatic 15 geothermal reservoir the overall seismic velocity structure mainly reflects the lithological heterogeneity of the system, whereas indicators for reservoir permeability and fluid content are deducible from the magnitude of seismic attenuation and the critical frequency at which the peak of attenuation and maximum velocity dispersion occur. The study demonstrates how numerical modeling provides a valuable tool to overcome interpretation ambiguity and to gain a better understanding of the hydrology of geothermal systems, which are embedded in a highly heterogeneous host medium.

Introduction
Magmatic geothermal reservoirs consist of permeable extrusive and intrusive rock formations, situated at depths where sufficiently high temperatures prevail. They are saturated with hot fluids, and usually heated by magma intrusions beneath the system. Evaluating the quality of such a reservoir requires an estimate of the fluid enthalpy and of the host rock permeability.
Seismic methods are amongst the most efficient exploration techniques to image the deep subsurface. The key quantities which 5 can be obtained from a seismic survey are the geometry of subsurface interfaces (e.g. lithological boundaries, faults, fracture zones), the P-and S-wavespeeds (V P and V S ) of various rock units, and the corresponding seismic attenuation characteristics. The latter is expressed by the inverse of the P-and S-wave specific quality factors Q P and Q S . The challenge in seismic interpretation is to link these seismic properties with the geological/hydrological properties of interest.
To constrain the seismic interpretation, it is recommended to measure the elastic and anelastic rock properties of small rock 10 specimens in the laboratory under in situ pressure-, temperature-, and fluid content-conditions. However, in magmatic geothermal systems, the host rock is often highly impermeable and the fluid transport predominately takes place within macro-fracture networks, rather than through the matrix. Such fractures are not present in the rock samples investigated in the laboratory, due to their limited size. Therefore, laboratory experiments only provide the properties of relatively intact rock and indicators for the presence or absence of fluids need to be deduced from fluid-rock interactions at larger scales through rock physics 15 concepts. Various such concepts of differing complexity have been used over the last 20 years to interpret seismic tomograms from geothermal exploration campaigns in magmatic environments.
Perhaps the simplest and most straightforward way to infer the presence of fluids in seismic interpretation is to recognize that V P is more sensitive to fluid saturation effects than V S , as the presence of liquids tends to increase V P but not significantly change V S . Thus it is common practice to deduce fluid saturation from seismic tomograms by interpreting the V P /V S ratio in 20 a qualitative manner. For instance Sanders et al. (1995) and Jousset et al. (2011) interpreted V P /V S anomalies to be indicative of the presence of supercritical fluids in a formation of the geothermal system in the Long Valley Caldera, California and in the Hengill volcanic complex in Iceland, respectively. Gunasekera et al. (2003), who conducted a time-lapse local earthquake tomography study in The Geysers, California, over a time period of 7 years, interpreted temporal variations in a V P /V S anomaly during the time of observation as an indication of water depletion resulting from reservoir operation. 25 For a more quantitative seismic interpretation, à priori information of the physical properties of mineral and fluid phases occurring at depth has to be taken into account. For instance, Julian et al. (1996) interpreted V P /V S anomalies observed in The Geysers, California, in terms of steam pressure, based on a mixing law of fluid and rock mineral properties. A more common way to incorporate fluid properties is through well-known fluid substitution theories, such as those of Gassmann (1951) and Biot (1956a and1956b), together with estimates of the rock frame mechanics (e.g., Nur and Simmons, 1969;Dvorkin et al., 30 1999). Husen et al. (2004) processed local earthquake tomography data from the Yellowstone volcanic field, Wyoming, while Vanorio et al. (2005) carried out a comparable study of data from Campi Flegrei, Italy. They concluded from fluid substitution calculations that V P /V S -ratio anomalies were caused by gaseous pore fluids. De Matteis et al. (2008) acquired tomograms in based seismic interpretation. Moreover, it needs to be recognized that seismic techniques cover a wide range of frequencies, from less than 1 Hz for local earthquake tomography, to more than 100 Hz for active seismic investigations, to the tens of kHzrange for sonic borehole tools, and up to 1MHz for piezo-electric pulse experiments in the laboratory. Thus, it is important to not only model the seismic response of fractured rock at the low-and high-frequency limit, but also at intermediate frequencies.
Different analytical approaches exist to account for velocity dispersion and attenuation due to wave induced fluid flow.

30
For instance Chapman (2003) describe the relaxation of fluid pressure between fluid inclusions, where the compliance of the inclusions is obtained from Eshelby (1957)'s theory. By contrast, Pride et al. (2004) and Gurevich et al. (2009) used Biot (1941's theory of poroelasticity, whereas Liu et al. (2009) used the theory of viscoelasticity to consider fluid flow in a double-porosity medium. Such theoretical models are based on some simplifying assumptions such as low fracture density or small elastic property contrasts, together with idealized geometries of heterogeneities. Motivated by this, numerical modeling approaches, based on the theory of poroelasticity as in the case of Masson and Pride (2007), Rubino et al. (2008), Wenzlau et al. (2010), and Quintal et al. (2011, became popular during the last decade, to complement analytical models. In this study, we use a numerical modeling technique, which is similar to those proposed by Rubino et al. (2008) andQuintal et al. (2011), to compute the seismic phase velocities and the frequency dependent wave attenuation in fluid saturated fractured 5 reservoirs. The reservoir is embedded in a magmatic-type environment, as it is typical for Iceland. We first define the physical properties of intact rocks based on the results of laboratory experiments reported in the literature. We take into account the diversity of typical rock types, which are shown to exhibit a large variability of hydro-mechanical properties. Then, for the up-scaling to the dimensions of macro-fractures, we study the properties of individual fractures in dependence of the hosting intact rock using a semi-analytical effective medium approach, which is based on Eshelby (1957)'s elastic field theory. Once the 10 parameters, which describe the physics of fractured rock volumes, are defined for ambient confining pressures under which the fractures are considered to be open, we study how each of these parameters depends on lithostatic stress, under which fractures close gradually. After this parameterization study, we finally apply the numerical model to a fractured geothermal reservoir in Iceland, as described in the structural geology literature. We examine how the frequency-dependent seismic properties of a rock containing a fracture network are affected by its saturating fluid, and how the observed fluid effects differ, depending on 15 the hosting lithology and on the effective lithostatic stress. To study the effects of fluids on the seismic properties of fractured rock, we use a numerical modeling technique which is based on the work of Quintal et al. (2011). It primarily involves Biot's (e.g. 1941) theory of poroelasticity and the principle of 20 conservation of linear momentum; where σ is the stress tensor, whose components in 2-D are related to the corresponding elements of the strain tensor by the constitutive law (2) 25 Here α = 1 − K d /K s is the Biot-Willis coefficient, λ = K d − 2/3G d is Lamé's constant, and δ ij is the Kronecker delta. The quantity K d is the drained frame bulk modulus, G d is the drained frame shear modulus, and K s is the bulk modulus of the solid (grain) phase of the porous rock. The drained state is equivalent to no fluid in the pores. The first two terms on the right hand side of (2) are consistent with Hooke's law of linear elasticity, whereas the additional term αP pore δ ij accounts for the stiffening of the rock in response to a pore pressure P pore . Biot (1941) completed his theory by adding the conservation of 30 fluid mass, under the assumption of fluid incompressibility. This requires that the flow rate into or out of an element of rock, described by Darcy's law for a global flow of liquid in a porous medium, is equal to the temporal change of fluid volume due to the deformation of the rock mass and due to the change of pore pressure. Transforming the mathematical formulation used by Quintal et al. (2011) into the space-frequency domain, the fluid transport equation is given by: where the imaginary quantity i and the angular frequency ω = 2πf represent the frequency domain-equivalent of the time derivatives. Quantity k is the hydraulic permeability, η f is the fluid viscosity, K f is the fluid incompressibility, and φ is the effective porosity.
To compute the poroelastic response of the medium, we simultaneously solve Eqs.
(1) to (3) for the stress relaxation resulting 10 from an imposed strain, using the COMSOL Multiphysics ® finite element solver. In its poroelastic representation on a finite element grid, a fractured rock as observed in nature, containing micro-fractures and macro-fractures of complex shape (Fig. 1a), is defined in a simplified manner (conceptual representation) by a composite of two poroelastic phases -the rock domain and the fracture domain (Fig. 1b). The rock domain represents the parts of the rock which are intact, apart from microscopic cracks which are not discretized individually, and it will be referred to hereafter as the intact rock. The fracture domain comprises 15 all the macroscopic fractures, which are in the model individually represented by smooth elliptic structures. We simply refer to them as fractures in what follows. As evident from (2) and (3), the hydro-mechanical behavior of each of these two media depends on a set of parameters, which are K d , G d , and K s , φ, and k for the solid phase of intact rock and fractures, and η f and K f for the saturating fluid phase. To distinguish between properties of the two media, we will mark intact rock properties with a hat superscript ('^') and the fracture properties with a tilde superscript ('~') throughout the text. 20 The model domain has undrained boundaries, meaning that there is no fluid flow across them. To conduct an oscillatory compressibility test, we simulate a vertical normal stress by a displacement disturbance ∆u in x 1 -direction to the top boundary, when referring to the coordinate frame in Fig. 2a, and we suppress any displacements in x 2 -direction at the left and right boundaries, and any displacement in x 1 -direction at the bottom boundary, as defined in Eq. (A1) in Appendix A. The stressstrain ratio resulting under these conditions yields the complex P-wave modulus, which is for a P-wave propagating towards For an oscillatory shear test, we apply a displacement ∆u in x 2 -direction to the top-boundary, and suppress any displacement in x 2 -direction at the bottom-boundary, while particles on the left and right boundaries are free to move in both directions x 1 and x 2 , as summarized in Appendix A by Eq. (A2). From the stress-strain relation calculated by this shear-test, we obtain the frequency-dependent complex shear-wave modulus for a S-wave propagating towards the x 2 -direction from the relation G c (ω) = 1 2 The angle brackets in Eqs. (4) and (5) denote the average over the entire modeling domain. Knowing the bulk density of the rock ρ b , seismic phase velocities can be obtained from the complex elastic moduli by (e.g., Casula and Carcione, 1992) and The attenuation factors are defined as the inverse P-and S-wave quality factors by (e.g., Casula and Carcione, 1992) The price for getting such a detailed characterization of seismic properties is that the method requires the determination 15 of a large set of parameters. In this study, we will give a detailed overview of typical values and likely ranges for each of these parameters, while accounting for the large diversity of rock types in magmatic geothermal systems. Furthermore, we will study how these values depend on lithostatic stress. A problematic feature with this approach is that parameterizing individual fractures by a homogeneous poroelastic medium and not by fluid filled cavities is a more conceptual rather than a direct physical representation. In particular, the definition of the fracture stiffness by intrinsic specific elastic moduliK d andG d neglects the 20 fact that in reality the elasticity of an open fracture is a complex interplay between the geometry of the void and the elasticity of the surrounding intact rock, which also involves a changing behavior of the intact rock due to the presence of the fracture, as has been described by Eshelby (1957).

Semi-analytical effective medium modeling
To study the properties of individual fractures under dry conditions, as required to determine the dry frame elastic moduliK d 25 andG d , we use the semi-analytical solution provided by the effective medium theory. The effective elasticity of a composite material, consisting of an isotropic elastic intact rock containing ellipsoidal inclusions which are filled with an isotropic elastic material (or empty as in our case), is calculated with the Mori-Tanaka method (Mori and Tanaka, 1973). An expression for the effective elastic tensor for the case where the ellipsoids are randomly distributed in a plane (holding one axis of the ellipsoids fixed, whereas the other two are randomly oriented as shown in Fig. 2b), is given by Pan and Weng (1995) as Here, C m is the elasticity tensor (in Voigt's matrix notation) of the intact rock, c f the volumetric concentration of inclusions, I the identity matrix and A is the eigenstrain concentration tensor, describing the strain under zero stress. The latter quantity is  Pan and Weng (1995) to be and In (11) and (12), angle brackets · denote the orientational average of the corresponding tensor, given in Appendix B.
Quantity C f is the elasticity tensor of the fracture filling-fluid phase and S is the Eshelby (1957) tensor, whose components for ellipsoidal inclusions can be calculated (e.g., Mura, 1987) from the aspect ratio of the ellipsoids and the elastic properties of the intact rock, i.e. fromK d andĜ d . Due to the random orientation of the ellipsoids in the x 1 -x 2 plane, the resulting effective elasticity tensor is transversely isotropic and the velocities of P-and S-waves propagating along the x 1 -axis are calculated from the corresponding components of the elasticity tensor C Eff by The effective medium theory is subject to several limitations in terms of the geometrical representation of fractured rock.
The underlying theory is exact only for non-interacting fractures (Kachanov, 1992), and is consistent with the upper and lower Voigt-Reuss bounds only in the case of low volumetric fracture density (Berryman and Berge, 1996). Furthermore, as stated by Kachanov (1992), the assumptions of non-interacting fractures and of small fracture density are not equivalent, since for non-20 randomly located fractures, the interaction might still be strong even for a dilute fracture density. For these reasons, fractures are considered to be well separated from each other and randomly located in space.

Dry fracture elasticity estimation
Compared with the benefits and drawbacks of the numerical modeling technique, the semi-analytical effective medium theory has complementary cons and pros. The effective medium is limited to non-interactive fractures, while the poroelastic theory 25 implemented on a finite element gird allows modeling the hydro-mechanical interaction of complex fracture networks. On the other hand, the effective medium theory implicitly includes the stiffness of the fractures, depending on the intact rock elasticity and the geometry of the fractures, whereas the numerical technique requires parameterizing individual fractures by a poroelastic medium, whereK d andG d are treated as fracture intrinsic material properties.
In the parameterization part of this paper, we will combine the two techniques to obtain appropriate values for the dry frame fracture stiffness'sK d andG d by varyingK d andG d until the stiffness of the overall fractured rock resulting from the poroelastic numerical modeling is consistent with that from the effective medium theory. To ensure that the 2-D numerical fractured rock model satisfies the requirements of the effective medium model, we generate models of randomly located, 5 randomly oriented and well-separated fractures of thin elliptic shape (black lines in Fig. 2a). The volumetric concentration of fractures is below 1%, what is below the limit for the low-fracture density assumption of 10-20% determined by Berryman and Berge (1996). We define an ellipse-shaped fracture-free area around each fracture (dotted ellipses in Fig. 2-a), and to guarantee that the fracture orientation is not biased by neighboring fractures, we successively place fractures within a circular area (dashed circle in Fig. 2a), which allows rotating the fracture by 360 • independently from the orientation of neighboring 10 fractures. As the numerical modeling is in 2-D, we assume the fractures to extend continuously in the out of plane direction over distances that are long compared to the in-plane fracture dimensions. Therefore, the 3-D effective medium model ( Fig.   2b) contains fractures with semi-major axis a 3 being much longer than semi-minor axes a 1 and a 2 , whereby the solution of the effective medium theory with ellipsoidal inclusions converts to one from a composite containing elliptic cylinders. The aspect ratio a 1 /a 2 , the fracture density c f , and the intact rock properties are chosen to match those of the numerical model. 15 When estimating values of the stiffnessK d andG d of dry fractures, no fluids are involved and the non-interaction condition is satisfied also for the numerical model in terms of fluid flow between fractures. Once appropriate values ofK d andG d are found, we will extend the complexity of the numerical fractured rock model beyond the capability of the effective medium theory, giving an example of modeling the seismic response of rocks containing fluid saturated interconnected fracture networks. Also for this case, the volumetric fracture density is still below 1% and, thus, the low fracture density-condition is still fulfilled. Here, 20 uncertainty arises from the fact, that the surrounding material of individual fractures also includes a small fraction of weaker material as fractures intersect. Uncertainties related to this effect can be reduced by using a more comprehensive effective medium theory than the one presented here, such as the self-consistent approach (Mavko et al., 2009, p. 185), which introduces a fractured background medium in an iterative manner.

25
In the present study, we focus on Icelandic geothermal systems. Iceland is a large subaerial part of the worldwide system of midocean ridges. Thus, the crust is to some degree of oceanic type, but anomalously thick with a maximum Moho depth of around 20 to 40 km (Darbyshire et al., 2000). The crustal sequence has been compared with the classical oceanic ophiolite sequence Bjarnason and Schmeling, 2009), but is of larger structural and chemical complexity (Gudmundsson, 1995). In a brief summary, the stratigraphy of the upper few kilometers of the Icelandic crust can be subdivided into four 30 lithological units (Fig. 3). At the shallowest depths, extrusive rocks dominate, forming interlayered sequences of (A) pyroclastic deposits (hyaloclastits, tuffs, scoria, etc.) and (B) lavaflow deposits (dense and vesicular basalts). In lower regions, (C) dyke and sheet intrusions (dolerite) become more and more abundant, which reach down to depths were (D) intra-crustal crystallized magma chambers (gabbro bodies) exist, which are found at depths as shallow as 1-2 km (Gudmundsson, 1995), but typically they occur at greater depth.
The physical properties of these rock types depend to some degree on their chemistry, which is predominantly of basaltic composition but, to a lesser extent, also magmatic rocks crystallized from intermediate and acid magmas exist (Gudmundsson, 1995). Depending on the temperatures and the intensity of fluid circulation through the formations, the chemistry of the rocks 5 is modified by hydrothermal alteration, what additionally increases the variety of minerals, each with potentially different physical properties. But not only the chemistry influences the physical properties of the rocks, also the rock texture has a strong influence. Dense gabbros are different from e.g. vesicular basalts or a highly porous tuff, independently from their chemical compositions. This also results in a large variability of the seismic properties as has been reported, e.g. by Vanorio et al. (2002) for pyroclastic rocks and by Grab et al. (2015) for basalts, dolerites and gabbros. Accordingly, a large variability is expected 10 for the properties of the intact rock in our models, and a large volume of data is required to provide well-grounded estimates.
Pyroclastic deposits are a typical feature of on-land volcanism. Lavaflow deposits, dykes and sheets, and magma chambers can also be found at submarine mid-ocean ridges. Therefore we can include the database of the ocean drilling programs for determining their physical properties.

Model parameterization for ambient lithostatic stress 15
The poroelastic model of fractured rock consists of two subdomains, intact rock and fractures. Their solid matrix are described by the same type of parameters. For the intact rock these are the dry frame bulk modulusK d , dry frame shear modulusĜ d , grain bulk modulusK s , dry bulk densityρ b , effective porosityφ, and hydraulic permeabilityk. The analogous parameters for the fractures areK d ,G d ,K s ,ρ b ,φ, andk. 20 We defined intact rock as those parts of the rock embedding the macroscopic fractures. Due to their limited size, rock samples investigated in the laboratory typically are free of such macro fractures, that is why we will refer to laboratory studies to determine intact rock properties. We here present a compilation of published results, which include more than 500 rock samples in total of diverse types from on-land volcanic systems as well as samples included in the database of the Deep Sea Drilling

Intact rock properties
Program and the Ocean Drilling Program.
25 Figure 4 illustrates in the form of cross plots values for all solid-constituent parameters as they have been reported in the literature. Assigning each rock sample to one of the main lithologies introduced in Section 3, we can study typical physical properties for each of these lithologies. This is depicted in Fig. 4a for the drained bulk modulus, with the boxes indicating the 25th and 75th percentiles, and dots are values outside these percentiles.
Values for the dry frame elastic moduli,K d andĜ d are obtained from velocities V P and V S , measured in the laboratory at 30 ultrasonic frequencies, from the relationŝ provided the wavespeeds were measured on dried rock specimens under drained conditions. In Fig. 4b,K d is plotted against G d , indicated by the open symbols, together with the elastic moduli of saturated rocks shown by the filled symbols. As for all other parameters, we seek to establish regression relationships using appropriate functions, to get representative values to 5 parameterize the fractured rock models. For the dry frame elastic moduli, the best-fit relationship was found to be: with p 1 = 1.4 × 10 5 , p 2 = 0.5117, and p 3 = −11.5 × 10 9 .
The grain bulk moduliK s were estimated using Gassmann's (1951) fluid substitution theory, which usesK s together witĥ K d to predict the bulk modulus of the saturated rock, whereasĜ d is assumed to be independent of liquid saturation conditions.

10
There is evidence for the validity of this supposition in the data shown in Fig. 4b, where the bulk moduli of saturated rocks (filled symbols) tend to be higher than those of dry rocks (empty symbols), and no significant increase is observed for the shear moduli.
Applying 4c. From the regression analysis, we find 20 with p 4 = 5.82 × 10 10 , p 5 = 3.86 × 10 −12 , p 6 = 8.22 × 10 08 , and p 7 = 3.99 × 10 −11 . To test the validity of these estimates, we useK s together with the effective porosityφ obtained from Eq. (18) introduced below, to predict the saturated bulk moduli from the dry bulk moduli by fluid substitution. The resulting saturated bulk moduli are indicated by the dashed line in Fig. 4b, which agrees well with the observed saturated bulk moduli given by the filled symbols, which includes numerous samples not being used for theK s -estimation. 25 Most researchers who measured seismic wavespeeds also documented the density and the porosity of the rock samples in their publications. Densities are plotted againstK d in Fig. 4d, and an exponential relationship is indicated, which yields from the regression analysis the best-fit function: with p 8 = 2628, p 9 = 1.72 × 10 −12 , p 10 = −2898, and p 11 = −1.38 × 10 −10 . Values for the effective porosities are shown 30 in Fig. 4e, and the regression analysis yields the best-fit relationship with p 12 = 0.85, p 13 = −1.13 × 10 −10 .
Since only a few authors measured the hydraulic permeabilityk and seismic wavespeeds together, we refer to different publications to estimate values fork. They are plotted against the bulk density in Fig. 4f and the best fit was obtained with 5 with p 14 = −11.80, p 15 = 4.26 × 10 −05 , p 16 = −3.60 × 10 −03 , and p 17 = 2.51 × 10 −03 .
Equations (15) to (19) fully describe the solid frame of the intact rock as a poroelastic medium. As evident from Fig. 4a, it covers a wide range of different lithologies. For the following modeling of the seismic properties of fractured rock, we select parameters for 4 characteristic models, each representing a different lithological unit which we will refer to as lithology A-D.
Lithology A represents a typical pyroclastic rock (blue dots in Fig. 4), lithology B a light lava flow deposit (red dots), lithology 10 C a typical dyke or sheet intrusive (yellow dots), and lithology D a relatively dense gabbro body (purple dots). Corresponding parameters for these four models are shown in Table 1.

Fracture properties
At ambient stress, fractures are assumed to be completely open, meaning that fracture walls are not in contact with each other and they can be represented by open ellipses (Fig. 1b). They are considered to be empty, i.e. containing no fault gauge, thus we 15 set the fracture porosity to a high value,φ = 90%. Furthermore, the mineral composition is assumed to be homogeneous across both the intact rock and the fracture subdomains, with the grain bulk moduli of the two subdomains being identical,K s =K s .
This also has the consequence that the mineral density is the same for both subdomains and the dry bulk density of the fracture Estimates for the dry frame elastic moduli of fractures,K d andG d , are obtained by testing what values ofK d andG d are 20 needed to obtain the same overall fractured rock stiffness's M Num and G Num from numerical modeling as the values M Eff and G Eff calculated using the effective medium theory. This test is conducted for a fractured rock model with well-separated non-interacting fractures, which allows comparing the results from the numerical modeling with those of the effective medium theory as discussed above in Section 2.3. As an example, M Num and G Num calculated for a model with intact rock properties corresponding to lithology B, and for fractures with an aspect ratio a 1 /a 2 = 400, are shown in Fig. 5a and b by the colored

30
From theory, it is expected that the bulk and shear moduli of dry fractures are of similar magnitude (e.g. Lubbe et al., 2008), Experimental results, however, indicate a ratio for dry fractures which is in fact small but larger than and Nakagawa (2013) found the ratio to lie in the range 1.7 ≤K d /G d ≤ 1.9. Therefore, we use here a small ratio which is larger than one,K d /G d ≈ 1.5, and we find under this constraint a pair ofK d andG d values which results in a good agreement (low RMS value) between the numerical modeling and effective medium result, shown by the red dot in Fig. 5a, b and c. This procedure is repeated for all lithologies A-D and for aspect ratios varying between 100 and 600, resulting inK d andG d values 5 shown in Fig. 5d and e and listed in Table 2.
The hydraulic permeability of open fractures is defined from well-established empirical relations reported in the hydromechanical literature. Based on calculations of laminar flow between two parallel walls, the volumetric flux through a fracture was described by the cubic law to scale with the cube of the aperture (Witherspoon et al., 1980;Zimmerman and Bodvarsson, 1996), leading to hydraulic permeability of the fracture defined as (Mavko et al., 2009) where the fracture aperture (width) is given by the hydraulic aperture e h . This was defined to be the aperture needed to explain the actually observed flow rate through a fracture with rough fracture walls in a parallel plate model. Thus, e h can be regarded as a parallel-wall equivalent aperture. Based on experimental observations, a formula for calculating e h was suggested by Barton et al. (1985) to be: where h is the average mechanical aperture of the fracture. JRC is the joint roughness coefficient with JRC = 2.5 indicating a very smooth fracture whereas a fracture with JRC = 20 is extremely rough (Barton and De Quadros, 1997). In our models, we use JRC = 15, assuming relatively rough fractures.

Model parameterization as a function of lithostatic stress 20
To study how the solid frame of the intact rock behaves with depth, we analyze their dependence on the effective confining pressure P , which is defined as the difference between the actual confining pressure (or lithostatic stress) and pore pressure P = P conf − P pore . For the individual fractures we consider the simplest case of an effective stress applied normal to the long fracture axis, given as the normal effective stress σ n .
5.1 Intact rock properties as a function of confining pressure 25 In the laboratory, V P and V S are usually measured at various confining pressures, to simulate the lithostatic stress conditions as a function of depth. Such datasets allow one to study the change of the bulk modulus and the shear modulus as a function of confining pressure, approximated by the second order Taylor series expansions, and (23) and (24), are shown in different colors in Fig. 6a and b, together with laboratory data indicated in gray. Values for the derivatives are given in Table 3.
Referring to experimental studies, hydraulic permeability of intact rock as a function of confining pressure has been described by a log-log relationship, e.g. by Lee and Farmer (1993), with ambient pressure permeabilityk 0 and the coefficient b indicating the curvature of the function or the slope when plotting log(k) versus log(P ).
To determine values of b for the four cases of lithologies A-D, we analyze the datasets which comprise permeability measurements for intact rock cores at varying confining pressures. They are shown in gray in Fig. 6c. For each dataset, a best-fit curve  Table 3. Resulting graphs fork(P ) are shown 25 color-coded in Fig. 6c.
The change in intact rock porosity resulting from a change of applied effective pressure was described by Jaeger et al. (2007).
They presented an expression for the change of porosity at a specific effective pressure P n due to an increment of effective pressure dP , which can be written as a function of the dry frame bulk modulus and the grain bulk modulus:

30
where the initial porosity and the dry frame bulk modulus are also functions of the effective confining pressure, given at the initial pressure P n−1 = P n − dP . Thus, the porosity at a given effective pressure P can be calculated stepwise using small pressure increments dP and updatingφ(P n−1 ) andK d (P n−1 ) at each step. The grain bulk modulusK s is assumed to be approximately constant at varying confining pressures. The dry bulk density of the intact rock varies according to the porosity variation,ρ b = (1 −φ(P ))ρ s , assuming the density of the mineral phase ρ s is 5 constant.

Fracture properties as a function of normal stress
The closure of natural unfilled fractures under normal stress was described by Bandis et al. (1983) as a function of specific normal and tangential fracture compliance. These quantities are related to the dry frame bulk and shear moduli by (Mavko et al., 2009), in the cases whereB n andB t are the compliances of dry fractures.
Referring to experimental observations, Bandis et al. (1983) described the fracture closure under normal stress being of hyperbolic form, becoming asymptotic to a small non-zero residual aperture. Based on their expressions, we calculate the fracture aperture as a function of normal stress by whereM d,0 is the dry frame P-wave modulus of the fracture and h 0 is the mechanical aperture, both at ambient normal stress as indicated with the zero subscript. The coefficient a is defined as the maximum fracture closure coefficient, being the factor relating zero stress aperture to the maximum aperture closure at very high stress, dh max = ah 0 , which we introduced to eliminate the specific compliance in the expressions of Bandis et al. (1983). Stress-dependent apertures resulting from (23) 20 are given in Fig. 7c for lithologies A-D, with the maximum and minimum values in each case which results from the aspect ratio-dependency ofM d,0 (Table 2).
An expression for the dry frame elastic moduli of the fractures can also be derived from the work of Bandis et al. (1983), leading to: 25 This equation gives the P-wave modulus as a function of normal stress, ambient-pressure elasticity, and fracture closure coefficient a. Thus,M d,0 is an intrinsic material property, which we can define without requiring any information about the absolute aperture such as h 0 or h(σ n ).
From Eq. (29), the effective bulk and shear moduli of the fracture can be calculated according to the relations: where ν is the Poisson's ratio, describing the ratio between bulk and shear modulus, We can either assume a constant Poisson's ratio ν = ν 0 , as observed at low normal stress, or use a stress dependent Poisson's ratio ν = ν(σ n ), in order to take stress-dependent effects into account, such as the higher resistance against shear, when rough fracture walls become interlocked during the closure at increased normal stress.
In the literature, the fracture compliance at varying normal stress has been experimentally investigated in terms of specific compliancesB n andB t and we cannot directly compare it with the outcome of Eq. (29), because of the scaling by the absolute 10 fracture aperture (Eq. 27). Instead, we can compare the stress-dependent fracture compliance normalized by the initial zero stress compliance, where the absolute fracture aperture cancels out: This is shown in Fig. 7a, where theB n /B n,0 ratios reported in the literature are displayed in gray. The normalized P-wave compliances for lithologies A-D, calculated according to the right-hand side of Eq. (32) and using a = 0.75, are shown in 15 color, where for each lithology two graphs are shown for the maximum and minimum values, depending on the aspect ratio.
Using a constant Poisson's ratio, normalized S-wave compliances are identical to the normalized P-wave compliances. They are given in Fig. 7b for Model A-D again color-coded and compared with the literatureB n /B n,0 ratios in gray. Using a fracture closure coefficient a = 0.75, we observe a similar behavior of the fracture compliance as reported in the literature, such as those of Nakagawa (2013) with its strong decrease in fracture compliance already at moderate values of σ n or those of Lubbe 20 et al. (2008) which only decrease relatively slowly as σ n is elevated to 60 MPa ( Fig. 7a and b).
At zero σ n , we considered the fractures being open and we used the cubic law (Eq. (21) for calculating the hydraulic permeability of fractures. At increased σ n , the fractures start to close and the two fracture walls come locally into contact with each other. Due to these contacts, as stated by Cook (1992), the reduction of permeability with increasing σ n is more rapid than the cube of the joint closure and the cubic law is not valid anymore. For this reason, Cook (1992) extended the cubic law, 25 yielding In Eq. (33), the first additional term leads to permeability reducing faster with increasing σ n than the cube of the fracture closure. The last term is the residual permeabilityk res , which incorporates the approximately constant permeability at very Assuming that the closure of the fracture is entirely compensated by a decrease of fracture void, whereas the area which is occupied by the material comprising the microscopic roughness remains constant, leads to a fracture porosity as a function of The grain bulk modulusK s is assumed to be approximately constant at varying confining pressures, as for the intact rock.
The dry bulk density of the fractures varies according to the porosity variation,ρ b = (1 −φ(σ n ))ρ s , assuming that the density of the mineral phase ρ s is constant.
6 Example -Fractured rock of variable depth and lithology 10

Model setup
After determining the hydro-mechanical properties of the intact rock and the fractures, the final task is to model the seismic properties of a rock mass containing an interconnected fracture network, which is saturated with a fluid of specific properties. Here we present a synthetic example using a model containing a fracture network which represents a highly fractured geothermal reservoir. The network geometry is based on the structural geology observations of Gudmundsson et al. (2002), 15 who examined a highly fractured palaeo-geothermal field associated with the Husavik-Flatey fault in northern Iceland. The network is embedded in a host rock consisting of basaltic lava flow piles, meaning that the petrography is similar to the one we presented above as lithology B. The original depth of the system is estimated to be approximately 1.5 km below the Earth's surface. Today, the overburden has been largely removed by erosion and the fracture network outcrops at the surface, preserved in the form of mineral filled veins. 20 Gudmundsson et al. (2002) described in detail the statistics of the network geometry in terms of the spatial fracture frequency, as well as the orientation, the width, and the length-to-width relationships of the fossil fractures. This gives a complete image of the fracture network as it is required to setup our model. This is done using a model generator, which places fractures randomly within the model domain, incorporating the fracture network statistics by weighting functions which are identical to the observations from the Husavik-Flatey fault. The resulting model is shown in Fig. 8 (a), together with distributions of the 25 fracture orientations (b), apertures (c), and a cross plot of fracture aperture a 2 versus fracture length a 1 (d). The gray line in Fig. 8d indicates an aspect ratio a 1 /a 2 = 400, which is the average aspect ratio observed by Gudmundsson et al. (2002).
The variability of fracture aperture and aspect ratios is not only considered in the model geometry but also when assigning the parameters of the poroelastic media representing the fractures. Fractures of larger aperture exhibit higher permeabilities according to Eq. (21), and fractures of larger a 1 /a 2 aspect-ratios are stiffer than those with small a 1 /a 2 -ratios as shown in 30 Fig. 5d and e. Corresponding ranges for the (normalized) fracture compliance are plotted against lithostatic stress in Fig. 7a and b, as they were computed from Eq. (29) using the ambient pressure stiffness's listed in Table 2 (A1) for the compressibility test, and those in Eq. (A2) for the shear test. Frequencies were varied over a wide spectrum from 10 −2 Hz to 10 6 Hz.

Results
Example results for the deduced seismic properties of the fractured rock are shown in Fig. 9, plotted as the P-wave modulus M (a), S-wave modulus G (b), inverse P-wave quality factors 1/Q P (c), and inverse S-wave quality factors 1/Q S (d) against the logarithmically scaled frequency. They were computed for lithology B, undergoing a lithostatic pressure of 15 MPa.
Comparing the elastic moduli with the corresponding attenuation graphs, we observe the typical behavior in accordance with 15 Kramers-Kronig dispersion relation (Mavko et al., 2009). At the same frequencies at which M and G are strongly dispersive with distinct inflection points, 1/Q P and 1/Q S reach their local maxima. In the example shown here, these frequencies, referred to hereafter as characteristic frequencies f c , are f c = 10 −1 Hz and the less prominent one f c = 10 4 Hz for the P-wave properties, and f c = 10 −1 Hz and most pronounced f c = 10 5 Hz for the shear wave properties, the latter having secondary peaks at around 10 1 Hz and 10 3 Hz. Norris (1993) linked the characteristic frequency with the diffusion length l d , over which 20 wave-induced fluid pressure diffusion takes place, by the relation where D is the hydraulic diffusion coefficient, D = k/η f (φ/K f + (α − φ)/K s ) −1 M/M sat , and where M/M sat is the ratio of the dry frame P-wave modulus and the undrained (saturated) P-wave modulus. The spatial scales of the fluid pressure diffusion during numerical oscillation tests can be best inferred from snap shots of the Darcy fluxes q, which are deduced from the local 25 permeability values k and the pore pressure gradients through the equation: Absolute amplitudes of the fluid fluxes ||q|| = q 2 1 + q 2 2 (with q i being the i-th component of the flux vector) occurring under compressional oscillations at frequencies of 10 −1 Hz and 10 4 Hz are shown in Fig. 10a and b, respectively. Absolute fluxes arising under shear oscillations at frequencies of 10 −1 Hz and 10 5 Hz are depicted in Fig. 11a and b, respectively.

30
At the low frequency of 10 −1 Hz, we observe increased fluxes inside all fractures (see zoom-plots in Fig. 10a and 11a), as well as within large areas of the surrounding intact rock. This shows that during one oscillation cycle, the pore fluid flows from the strongly compressed compliant fractures deeply into the stiffer, and less permeable, intact rock. Thus it is a flow between heterogeneities, with the stiffness and permeability values differing by several orders of magnitudes. It takes place at scales larger than the pore scale but smaller than the wavelength, why this dispersion mechanism is commonly called the mesoscopic flow (MF) mechanism (e.g. Müller et al., 2010). For the example shown here, fluxes are more widespread for the compression experiment than for the shear experiment, what is in agreement with the larger 1/Q P -magnitude compared with the magnitude of 1/Q S at 10 −1 Hz, shown in Fig. 9c and d, respectively. Furthermore, the attenuation peak due to MF occurs 5 as a single maximum without minor side peaks. This is because the characteristic frequency of MF is predominantly controlled by the medium with the lower fluid mobility k/η f (Quintal et al., 2014), which is the intact-rock subdomain here, having a constant permeability throughout the entire model domain of 7 × 10 −18 m 2 . Solving Eq. (35) for the diffusion length using the poroelastic parameters of the intact rock and f c = 10 −1 Hz, we find l d ≈ 0.02 m, which is in good agreement with the width of regions with increased fluid fluxes in Figs. 10a and 11a.

10
At increased frequencies, fluid fluxes into the intact rock become less pronounced (Figs. 10b and 11b), because the shorter oscillation cycles limit the pressure relaxation by fluid flow from the fractures into the intact rock with its low permeability. In the intact rock, fluid flow only takes place within direct vicinity of the fractures and is more pronounced at the tips of individual fractures. Thus, as frequency increases, fluid flow concentrates more and more in the highly conductive fractures. This flow is driven by pressure gradients between different interconnected fractures, which undergo various degrees of compression, either 15 because they are oriented differently relative to the direction of the applied oscillation stress, or because they are of different stiffness. As we identify fluid flow between different fractures at these higher frequencies, the corresponding dispersion mechanism is equivalent to the squirt flow (SF)-type (e.g. Müller et al., 2010), but at a larger spatial scale. The stress for the compressibility test was oriented along the x 1 -axis, when referring to the coordinate frame in Fig. 8a, and the stress of the shear experiment was parallel to the x 2 -axis and being of the dextral form. Therefore, fluid flow dominates in different fractures for 20 the compressibility and for the shear test, which is why there are different numbers of peaks in the 1/Q P -and 1/Q S -plots, and why they are occurring at different frequencies varying between around 10 −1 and 10 5 Hz. This range can be explained by the fact that the characteristic frequencies linearly scale with fracture permeability according to Eq. (35), which in the case shown here are within the range of 10 −8 and 10 −11 m 2 . The diffusion length required to explain the observed characteristic frequencies is l d ≈ 0.1 m which is similar to the half-length of most fracture segments between fracture-intersection points in 25 our model (Fig. 8a).
Next, the seismic properties of all lithologies were modeled. The numerical results obtained for lithostatic pressures of 15 MPa are shown in Fig. 12. The red graphs, representing lithology B, are identical with those shown in Fig. 9, whose characteristics were discussed above. Comparing first the absolute magnitudes of the elastic moduli of all four lithologies, we observe that the fractured rock models A to D become successively stiffer, as a logical consequence of the increased 30 stiffness of all the rock components. The attenuation peak, which we interpreted to be of MF-type, occurs at characteristic frequencies f c ≤ 10 1 Hz. Comparing this peak for the four lithologies A-D, we observe a decrease in amplitude, A A > A B > A C > A D , which can be explained by the intact rock porosities being strongly decreasing in the orderφ A >φ B >φ C >φ D .
Higher porosities entail larger amounts of fluid in the saturated rocks and, hence, more energy is consumed by fluid flow leading to higher attenuation. This effect opposes and dominates over the effect of varying stiffness contrasts, which here are of similar magnitudes for the four lithologies, but which would amplify the attenuation when the stiffness contrasts increase. Regarding the characteristic frequency, it is observed to decrease for the four models, f c,A > f c,B > f c,C > f c,D , which is in agreement with Eq. (35) and the fact that the intact rock hydraulic permeability of the four lithologies progressively decreases ask A >k B >k C >k D . The opposite effect is observed for the attenuation peaks at the higher frequencies, which can be related 5 with the SF-type mechanism. For the lithologies A-D the peaks occur with increasing amplitudes A A < A B < A C < A D and at increasing characteristic frequencies f c,A < f c,B < f c,C < f c,D . This is because there is less resistance against the fractures closure under lithostatic stress for the softest rock of type A compared to B, C and subsequently D (see Fig. 7c). Therefore, fractures embedded in lithology D remain the most open and retain their original permeability and fluid saturated pore space the most, which is why the SF-attenuation peak is largest and occurs at the highest frequency for lithology D and subsequently 10 lowers for lithology C, B, and A.
Finally, to study the effect of lithostatic stress on the seismic properties of fractured rock, the elastic moduli and the seismic attenuation are computed for lithology B, on which an effective lithostatic stress is applied ranging from ambient pressures up to 120 MPa. Numerical results are shown in Fig. 13. We observe the elastic moduli M and G to strongly increase with stress, which is predominantly due to the strong stress-dependence of the fracture stiffness ( Fig. 7a and b), together with the 15 less dominant stiffening of the matrix ( Fig. 6a and b). Comparing the attenuation peaks for varying stress, we observe the MF-peaks occurring at an approximately constant frequency of 10 −1 Hz, which is consistent with the approximately constant permeability of the intact rock, varying by not more than one order of magnitude for a confining pressure ranging from 0 to 120 MPa. The SF-attenuation peak, however, occurs at a characteristic frequency of f c > 10 6 Hz in the case of zero lithostatic stress, and decreases down to f c ≈ 10 3 Hz for a lithostatic stress of 30 MPa. For higher lithostatic stress the characteristic 20 frequency decreases even more, down to frequencies which overlap with those of the MF-type dispersion, where the peaks start to overlap indistinguishably. For both types of attenuation mechanism, the amplitudes decrease with increasing lithostatic stress. This is due to the reduced porosity, and, hence, a reduced pore water content, combined with a stiffening of both, intact rock and fractures at elevated lithostatic stress. This gives rise to lower amplitudes of the wave-induced fluid pressure gradients due to the smaller compressibility contrast. As a consequence of these two effects, the amount of fluid flow is reduced and less 25 energy is dissipated, leading to smaller attenuation peaks at increasing lithostatic stress.

Discussion and Conclusion
In fluid-saturated fractured rock, viscoelastic interaction between the intact rock, the fractures, and the saturating pore fluid causes velocity dispersion and seismic wave attenuation. The underlying mechanisms have been studied in the past by various researchers, as summarized by Müller et al. (2010), and there is a broad consensus about how the degree of seismic wave 30 attenuation and the characteristic frequency at which it occurs depends on the hydro-mechanical properties of the materials constituting the rock. Petrophysical models which consider such viscous fluid flow are able to link seismic quantities which are measured in geothermal exploration campaigns with the hydrological properties. The reason why these models have not been used routinely to date in seismic interpretation is to a large extent because they depend on many input parameters, some of which are difficult to quantify.
We have determined the input parameters for magmatic geothermal systems, as required in numerical oscillation tests, and wide ranges were observed for most properties when considering the high diversity of magmatic rock types. Most of the input parameters also depend on lithostatic stress, which is why we provided a compilation of functions to calculate the input 5 parameters for varying effective stress. Using these parameters, we computed the seismic properties of rock volumes containing an interconnected fracture network saturated with liquid water. Results from the numerical modeling demonstrate how seismic velocities and attenuation factors strongly depend on the lithology. This was already established for P-and S-wave velocities in our earlier experimental study (Grab et al., 2015). Here, this is ground-truthed by a large database extracted from the literature, which shows that the seismic velocity structure of magmatic geothermal systems primarily reflects the subsurface lithology.

10
The effects of reservoir permeability and fluid content are only minor. Interpreting seismic data in terms of hydrological target parameters against the contrast of this background heterogeneity can be achieved by studying the seismic attenuation, i.e., the decrease of seismic amplitudes with increasing distance of travel, and if available, by studying the velocity dispersion, viz., seismic velocity differences at varying frequencies.
Our modeling results show how the magnitudes of seismic attenuation and its dispersion are associated with stiffness con-15 trasts and porosity. Large attenuation peaks were found for a rock volume containing a network of open fractures, which decreased considerably when subjecting the fracture network to elevated lithostatic pressures forcing the fractures to close.
The characteristic frequency, at which the attenuation reaches its peak, is linked with the fluid mobility, which is a measure of hydraulic permeability and fluid viscosity. At low seismic frequencies, the attenuation is observed to be controlled by mesoscopic fluid flow from fractures into the surrounding porous intact rock, with the characteristic frequency linearly scaling with 20 intact rock permeability and fluid viscosity. At sonic up to ultrasonic frequencies, attenuation is associated with squirt flow between interconnected mesoscopic fractures which are compressed to differing degrees during normal and shear oscillations.
Here, the characteristic frequency linearly scales with fracture permeability and fluid viscosity.
The spread in the observed critical frequencies illustrates that fluid effects in fractured rock can be detected with various seismic techniques (passive, active, sonic, etc.) or in the ideal case by the combined use of different seismic techniques to cover 25 a broader frequency spectrum. On the other hand, there seems to be no general rule that governs the frequencies at which the conditions are given to assume either the relaxed state (low frequency limit) or unrelaxed state (high frequency limit), beyond which traditional rock-physics concepts are strictly valid. Thus, concepts which incorporate wave-induced fluid flow, like the one we presented in our study, can help improve the quantitative interpretation of all kinds of seismic data.
In the scope of this study, we modeled the influence of wave-induced fluid flow on seismic properties for the wide range Experimental investigations on the stiffness and permeability of intact rock and fractures at elevated temperatures are rare. In 5 general, it is known that an increase in temperature results in softening of the intact rock. Thus, increasing the temperature may cause similar behavior to moving from a stiff lithology to a more compliant one as examined in this study.

Appendix A: Boundary conditions
The boundaries Γ of the model domain Ω consist of undrained boundaries. To conduct an oscillatory compressibility test, we simulate a normal stress by a displacement disturbance ∆u in the x 1 -direction to the top boundary Γ T , when referring to the 10 coordinate frame in Fig. 2a, and we suppress any displacements in the x 2 -direction at the left Γ L and right Γ R boundaries, and any displacement towards the x 1 -direction at the bottom boundary Γ B , i.e. rigid boundaries at right, left, and bottom, given by Accordingly, we apply a displacement in x 2 -direction to Γ T for the oscillatory shear test, and suppress any displacement 15 towards the x 2 -direction at Γ B , Meanwhile, particles on Γ T and Γ R are free to move into both directions x 1 and x 2 .

Appendix B: Orientational average
The orientational average of a fourth-order tensor is where the rotation around the third principal axis x 3 by the angle θ is obtained by applying the transformation law with R being the corresponding entry of the rotation matrix This leads to the following expressions for the averaged elasticity tensor X , given in Voigt's matrix notation: X 11 = X 22 = 1 8 (3X 11 + 3X 22 + X 12 + X 21 + nX 66 ) where the factor n depends on the definition of the shear components X 44 , X 55 , and X 66 , when transforming the fourth rank 5 elasticity tensor into Voigt's matrix notation.