the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
Seismic attenuation and dispersion in poroelastic media with fractures of variable aperture distributions
Nicolás D. Barbosa
J. Germán Rubino
Beatriz Quintal
Considering poroelastic media containing periodically distributed parallel fractures, we numerically quantify the effects that fractures with variable aperture distributions have on seismic wave attenuation and velocity dispersion due to fluid pressure diffusion (FPD). To achieve this, realistic models of fractures are generated with a stratified percolation algorithm which provides statistical control over geometrical fracture properties such as density and distribution of contact areas. The results are sensitive to both geometrical properties, showing that an increase in the density of contact areas as well as a decrease in their correlation length reduce the effective seismic attenuation and the corresponding velocity dispersion. Moreover, we demonstrate that if equivalent physical properties accounting for the effects of contact areas are employed, simple planar fractures can be used to emulate the seismic response of fractures with realistic aperture distributions. The excellent agreement between their seismic responses was verified for all wave incidence angles and wave modes.
 Article
(3725 KB)  Fulltext XML

Supplement
(5919 KB)  BibTeX
 EndNote
Fractures in rocks occur in a wide range of scales (from microscale to continental), and their identification and characterisation are important tasks for several areas such as oil and gas exploration and extraction, production of geothermal energy, nuclear waste disposal and civil engineering works, among others (Schoenberg and Sayers, 1995; Metz et al., 2005; Tester et al., 2007). Given that seismic waves are known to be significantly affected by the presence of fractures (e.g. anisotropy, attenuation, dispersion, scattering), seismic methods are a valuable tool for detecting and characterising fractures. An important cause of seismic attenuation and velocity dispersion occurs when a fluidsaturated heterogenous rock is deformed by a propagating wave. In such a case, the compressibility contrast between the heterogeneities creates fluid pressure gradients. Then, the fluid pressure returns to equilibrium through a fluid pressure diffusion (FPD) process during which energy is dissipated due to viscous friction (e.g. Pride and Berryman, 2003; Gurevich et al., 2009; Müller et al., 2010). In the case that the rock heterogeneity is a fracture, FPD is generated as a consequence of the compressibility contrast between the fracture and the embedding background. Furthermore, when a propagating wave deforms a medium having intersecting and fluidsaturated fractures, a pressure gradient can also arise between them, resulting in additional energy dissipation due to FPD but affecting higher wave frequencies than the previously mentioned phenomenon (Rubino et al., 2013; Quintal et al., 2014, 2016).
Several authors have studied fracturerelated FPD effects on seismic attenuation and velocity dispersion (Chapman, 2003; Brajanovski et al., 2005; Carcione et al., 2013; Rubino et al., 2013; Quintal et al., 2014; Caspari et al., 2016, 2019), as well as on the effective anisotropy (Maultzsch et al., 2003; Masson and Pride, 2014; Tillotson et al., 2014; Amalokwu et al., 2015; Rubino et al., 2016, 2017) and scattering (Nakagawa and Schoenberg, 2007; Barbosa et al., 2016; Caspari et al., 2019), based on experimental or theoretical works. A common approach to study seismic attenuation and velocity dispersion in fluidsaturated fractured media consists of numerically solving Biot's (Biot, 1941, 1962) poroelasticity equations (Masson and Pride, 2007; Quintal et al., 2011, 2014; Rubino et al., 2013). In such studies, fractures are modelled as very compliant, highly porous and highly permeable heterogeneities embedded in a much stiffer, less porous and less permeable homogenous background (Pride et al., 2004).
Fractures can be conceptualised as two uneven surfaces in contact, which produce variable separation between their boundaries or walls (e.g. Montemagno and PyrakNolte, 1999; Jaeger et al., 2007; Masson and Pride, 2015). In spite of this, fractures are often modelled as simple thin layers of constant thickness. The validity of such representation relies on assumptions such as incident wavelengths larger than the fracture microstructure, an appropriate selection of the properties of the material filling the fractures, among others. Furthermore, analytical solutions are based on simplifications such as assuming elastic materials, idealised geometries for the fracture microstructures and/or neglecting heterogeneities interactions. In this sense, Hudson and Liu (1999) obtained equivalent mechanical properties for a constantthickness layer by representing a fracture as a plane distribution of circular dry or fluidfilled cracks (i.e. open regions of fractures). They found that the geometrical properties that dominate the mechanical behaviour of the fracture are cracks density and length. Still, their model assumes elastic media and crack or contact areas interactions are neglected. Also in an elastic framework, Zhao et al. (2016) have demonstrated the important role played by cracks' interaction on the stress field and, consequently, their effects on the overall rock stiffness. Hudson et al. (1996) proposed an analytical solution for obtaining the equivalent elastic properties of a constantthickness layer by representing a fracture as an aligned distribution of circular cracks, accounting for their interaction and considering fluid or viscoelastic material filling the cracks. However, this model is limited to an elastic background and a single fracture. By considering 2D models and following a poroelastic approach based on Biot's (Biot, 1941) equations, Rubino et al. (2014) numerically studied seismic attenuation and dispersion in fluidsaturated fractured media considering several distributions and densities of fracture contact areas. They found that contact areas have a strong effect on the level of seismic wave attenuation and dispersion caused by FPD between fracture and background, and also that the distribution of contact areas and material properties plays an important role on the effective seismic response. Nevertheless, they performed 2D simulations which cannot account for realistic distributions of fracture aperture and contact areas. To our knowledge, the analysis of the impact of realistic aperture distribution of fractures on seismic attenuation and dispersion due to FPD remains largely unexplored.
Nolte and PyrakNolte (1991) presented a stratified percolation workflow which can be used for generating fractures with realistic aperture distribution. The proposed methodology allows to define the structure of the fractures in terms of density and distribution of contact areas. PyrakNolte and Morris (2000) studied the effects that such geometrical fracture properties have on the fracture stiffness and on the fluid flow through the fractures under normal stress changes. Moreover, they established a relation between the hydraulic and mechanical behaviours of fractures. However, FPD effects on the seismic response of fractures with realistic aperture distribution were not considered in their work.
In this work, we follow the workflow proposed by Nolte and PyrakNolte (1991) to generate fracture models and use them to numerically quantify the effects that fractures with variable aperture distributions have on seismic wave attenuation and stiffness modulus dispersion for a medium containing a periodic distribution of fractures. For normally incident compressional waves, we perform a sensitivity analysis in terms of density and distribution of contact areas. Subsequently, we extend the results for all incidence angles and wave modes by computing the stiffness matrix coefficients. Finally, we numerically show that the same frequencydependent seismic attenuation and velocity dispersion can be obtained for fractures with intricate aperture distributions and for fractures with constant aperture, provided that appropriate fracture properties are used in the latter case. These equivalent fracture properties must take into account not only the properties of the fracturefilling material but also the mechanical effects associated with the aperture distribution.
To study attenuation and dispersion of seismic waves in a fluidsaturated rock with parallel and periodically distributed fractures, we model fractures as poroelastic media embedded in a homogenous poroelastic background. The aperture of the modelled fractures can be spatially variable. In the present work, we refer as open regions of the fracture to the zones where the fracture walls are not in contact (nonzero aperture) and are filled with a highly permeable and porous material. Contact areas (zero aperture), on the other hand, are represented by a porous material having the same properties as the background medium. We define the density of the contact areas as the ratio between the area of the fracture walls in contact and the area of the entire fracture. In this work, we use similar material properties to those employed by Rubino et al. (2014) but assuming a lower permeability for the background (Table 1). The material properties of the background are representative of sandstone (Bourbie et al., 1987).
Assuming that the prevailing wavelengths are much larger than the fracture aperture and spacing, we can obtain the effective seismic properties of the fractured medium by performing oscillatory relaxation tests on a representative elementary volume (REV) of the medium. We solve the quasistatic poroelastic equations given by Biot (1941) with a finite element scheme. The equations are written in the up (denoting solid displacement and fluid pressure, respectively) formulation, as proposed by Quintal et al. (2011), but in the space–frequency domain following the approach of Rubino et al. (2009). The effective stiffness moduli are obtained by applying homogeneous harmonic displacements normal to a boundary of the model in the case of P waves, or a shear displacement in the case of S waves. The numerical model, which consists of one horizontal fracture embedded in a homogenous background (Fig. 1a, left), is the REV of a medium containing parallel and periodically distributed fractures. On the right side of Fig. 1a, we provide a sketch of the described test for the case of normal incidence of a P wave, in which a vertical displacement is applied to the top boundary while no normal solid displacements are allowed at the lateral boundaries and at the bottom of the model. Additionally, the model is fully saturated with water and the test is performed at undrained conditions; that is, no fluid flow across the model boundaries is permitted.
2.1 Equations of poroelasticity
Biot's (1941) equations in the space–frequency domain and in absence of external forces are
where σ is the total stress tensor, w is the relative fluid displacement vector, p is the fluid pressure, ω is the angular frequency, and i is the complex unity. The material properties κ and η are the permeability and the fluid viscosity, respectively. Equations (1) and (2), the total balance of forces and Darcy's law, respectively, are also known as consolidation equations. These equations are coupled through the constitutive equations of the porous medium (Biot, 1962):
where ϵ is the strain tensor, I denotes the identity tensor, and e is the cubical dilatation given by the trace of the strain tensor. The symbols λ_{m} and μ_{m} denote Lamé's parameters of the dry frame, and α and M are Biot's poroelastic parameters given by
where K_{m}, K_{s} and K_{f} are the bulk modulus of the dry frame, the solid grains composing the rock frame and the fluid, respectively, and (ϕ) is the porosity. Combining Eqs. (2) and (4), the following expression can be obtained:
where u is the solid displacement vector. Finally, Eqs. (1) and (7), with the total stress tensor given by Eq. (3), can be used to express the consolidation equations in terms of the unknowns u and p. In 3D, the solid displacement vector u has three directional components and the fluid pressure is a scalar.
2.2 Effective properties
We solve the system of Eqs. (1) and (7) employing Eq. (3) and using the finite element method. Considering auxiliary functions to represent the variables between nodes, the “weak” formulation is obtained by combining the differential equations (“strong” form) with natural (undrained) boundary conditions in an integral form, which allows for reducing the order of the spatial derivatives (Quintal et al., 2011). We input directly the weak formulation in the finite element software, COMSOL Multiphysics, where the boundary conditions corresponding to the described relaxation tests are imposed. The numerical model is discretised in tetrahedral elements (Fig. 1b), where quadratic interpolation is applied. The unknowns u and p are calculated for each considered frequency, from which the stress and strain fields are obtained in each element of the numerical model. For the case of a compressional wave propagating normal to the fractures (z direction), the complex effective Pwave modulus (H) and the Pwave attenuation, estimated as the inverse of the quality factor (O'Connell and Budiansky, 1978), are computed according to the following expressions:
where 〈σ_{zz}(ω)〉 and 〈ϵ_{zz}(ω)〉 represent the volumetric averages of σ_{zz} and ϵ_{zz} over the whole numerical model for each frequency (e.g. Masson and Pride, 2007; Rubino et al., 2009; Jänicke et al., 2014), and Re and Im correspond to the real and imaginary parts of a complex number.
In order to numerically analyse the general impact that fractures with variable aperture produce on seismic attenuation and velocity dispersion, we first consider fractures with simple geometries and distributions of contact areas. Then, we extend the investigation to fractures having realistic aperture distributions and perform a sensitivity analysis of the effective seismic response of the fractured medium in terms of density and correlation length of contact areas.
3.1 Fractures with simple aperture distributions
We first consider simple fracture models for illustrating general effects of contact areas distributions on the Pwave modulus normal to the fractures and the associated seismic attenuation. The numerical model is a cube of 4 cm sides having a horizontal fracture crossing its centre, that is, normal to the vertical (z) direction. This model is an REV of a medium containing a periodic distribution of parallel fractures with 4 cm separation between the central planes of two consecutive fractures. Figure 2 shows a representation of the central plane of a fracture with regular distribution of contact areas and another one with pseudorandom distribution. The aperture of the open regions of the fracture is constant and equal to 0.4 mm (blue regions), whereas the white square zones correspond to contact areas (i.e. aperture equal to zero). In both cases, the contact area density is 20 %.
The numerical results for the real part of the effective Pwave modulus and seismic attenuation normal to the fracture are presented in Fig. 3. When the P wave compresses the fractured medium, it creates a fluid pressure gradient between the fracture and the background due to their compressibility contrast. Consequently, a FPD process tends to equilibrate the pressures resulting in energy dissipation because of viscous friction in the pore fluid. Both models exhibit significant Pwave modulus dispersion and attenuation due to FPD between the fracture and the background. Although the responses presented in Fig. 3 are similar, the highest dispersion and attenuation is observed for the pseudorandom distribution of contact areas. The largest difference in their Pwave modulus responses occurs at the lowfrequency limit. We compare these results with the analytical solution of White et al. (1975) for a fracture represented as a thin layer of constant thickness filled with the same soft material but without contact areas. Two models having fractures of constant thickness are considered, one of a thin layer having constant aperture of 0.4 mm and another one having equal total volume of open regions (that is, removing the volume occupied by the contact areas) as in the models shown in Fig. 2. We observe that the presence of contact areas reduces seismic attenuation as they increase the fracture stiffness. For a thin layer with the same volume as the open regions of fractures of Fig. 2 (i.e. a constant thickness of 0.32 mm), we observe minor changes in the seismic response with respect to the thin layer having 0.4 mm aperture. This means that the effects of contact areas on the mechanical behaviour of the fracture are more significant than the effects of reducing the volume of the open regions.
For better understanding the impact of contact area distributions on the seismic responses shown in Fig. 3, we plot the real part of the vertical component of the total stress field at the low and highfrequency limits at the centre of the fracture (Fig. 4). The interaction between contact areas in the pseudorandom case results in a stress shielding, which means that the proximity of the contact areas caused a reduction of the stresses on them (Zhao et al., 2016). From Eq. (8), and given that the overall strain in the sample is the same in both cases, it follows that a reduction in the overall stress of the sample translates into a decrease of the effective Pwave modulus, in comparison with the regular distribution. At the lowfrequency limit, as there is enough time for FPD between the fracture and the background, the stiffening effect of the fluid in the fractures is minimal. Therefore, given that the magnitude of the interaction effects depends on the compressibility contrast between contact areas and regions of open fracture, these are maximal at low frequencies. Further, Fig. 3 shows that the real part of the Pwave modulus for both models converges to a similar value at the highfrequency limit. This occurs because there is no time for fluid pressure exchange between the fracture and the background during a halfwave period, and therefore the stiffening effect of the fluid saturating the fracture is maximal. Hence, the effects of the interaction of the contact areas on the stress field are minimal, since the compressibility contrast between the background and the fracture is reduced (Fig. 4).
3.2 Fractures with realistic aperture distributions
In order to analyse the seismic response of realistic fractures, we perform numerical simulations considering fractures with variable aperture distributions generated following the stratified percolation approach of Nolte and PyrakNolte (1991). Using this approach, a realistic spatial distribution of fracture contact areas can be generated with controlled statistical properties such as correlation length and density. Moreover, it allows us to produce variable aperture distributions of the open regions. A description of how the fracture models are generated is given in Appendix A. The consideration of such aperture distributions as realistic is based on their comparison with the imaging of the aperture distribution of a natural fracture network presented by Montemagno and PyrakNolte (1999). The aperture of the fractures is given by the distance between the walls which are perfectly symmetrical with respect to the centre of the fracture (Fig. 1). After generation of fractures, the correlation length of contact areas is calculated following the approach of Blair et al. (1993), also used by PyrakNolte and Morris (2000), which represents an approximation to their mean length. Consequently, for a given contact area density, as the correlation length decreases, the fracture exhibits more contact areas with smaller sizes and a narrower distribution of distances between them. Thus, increasing the correlation length of contact areas produces an increase in the mean distance between contact areas, that is, the mean length of the open regions. Using these kinds of fracture models, PyrakNolte and Morris (2000) discussed the effects that contact area distributions produce on the mechanical properties of a fracture (i.e. specific stiffness), and they showed that uncorrelated distributions of contact areas produce stiffer fractures than correlated ones. This is in agreement with the results presented in Fig. 3, since a regular distribution of contact areas is analogous to an uncorrelated distribution considering that both have a narrow distribution of distances between contact areas.
Figure 5 shows four fracture aperture distributions chosen for analysing the effects of density and correlation length of contact area on the seismic response of a fractured medium. To do so, we consider four REVs as the one presented in Fig. 1 with a different distribution of apertures for the horizontal fracture given by the models A, B, C and D of Fig. 5. The real part of the Pwave modulus and attenuation for a normally incident wave are shown in Fig. 6. We observe that a higher contact area density as well as a lower correlation length (model C) produce a stiffer fracture, as it can be seen from the relatively high and nondispersive real part of the Pwave modulus. These observations can be understood by looking at Fig. 4 as, in this case, a regular distribution of contact areas is equivalent to an uncorrelated distribution. Moreover, these results are in agreement with analytical solutions, considering that an increase in the correlation length represents also an increase in the mean crack (or open region) length (Hudson and Liu, 1999; Zimmerman and Main, 2004). In such solutions, the excess compliance of a fracture increases with a larger crack length and decreases with a greater contact area density. Fracture B has the lowest density and the highest correlation length of contact areas. As expected, it is effectively the most compliant as can be seen in Fig. 6, given that it exhibits the lowest Pwave modulus and highest dispersion and attenuation. Interestingly, fractures A and D show comparable seismic responses despite their different aperture distributions. This is related to the fact that increasing the correlation length compensates for the increase in contact area density, resulting in fractures with similar effective compliances. Lastly, note that although the effect of the distribution of contact areas is maximal at the lowfrequency limit, it also plays an important role for the effective compliance of the rock at the highfrequency limit (Fig. 6).
A common assumption in analytical models is that fracture compliance depends on the fracture volume and the distribution of the fracture microstructure, such as mean crack radius (e.g. Guo et al., 2017). For analysing the effect of the variable aperture in the open regions of the fractures, we consider simpler fracture models having the same contact area distributions of those shown in Fig. 5 but setting constant the aperture in the open fracture regions. The aperture is fixed to their mean value (i.e. 0.4 mm) (Fig. 7). We refer to this process as binarisation of the fractures, since it results in two values for the fracture aperture: zero in the contact areas and 0.4 mm in the open regions. During this process, the volume of the open regions, that is, the volume of the most compliant poroelastic material, remains unchanged. The seismic responses are also shown in Fig. 6, and we observe good agreement between results from fractures with binarised aperture (dashed lines) and those from the previous fractures with variable aperture in the open regions (solid lines). From this analysis of Fig. 6, the volume of the open regions of the fractures and the distribution of contact areas are the main characteristics controlling the fracture seismic response, which is in agreement with analytical models.
In the analysis presented above, each fracture is one realisation of a pseudorandom generation process with given contact area density and correlation length. In order to analyse the variability of the results, we generated several fracture models with the same characteristics as model B (Fig. 5), which are referred to as realisations. We choose model B because it exhibits the highest attenuation (Fig. 6), and hence the variability of the results is expected to be higher than for uncorrelated distributions of contact areas (such as models A and C). For numerical convenience and supported by the comparison shown in Fig. 6, we consider binarised fracture models as in Fig. 7. The fracture models illustrated in Fig. 8a have equal contact area density (5 %) and correlation length (2.8 mm) to fracture B. In Fig. 8b, the corresponding real part of the Pwave modulus and attenuation are displayed as black circles for frequencies representative of the relaxed and unrelaxed limits (1 and 100 kHz, respectively) and at the attenuation peak frequency (∼100 Hz). The seismic response for model B of Fig. 5 is plotted as a solid red line for reference. The standard deviations of the real part of the Pwave modulus normal to the fractures are presented in Fig. 8c as a function of the number of realisations. We observe larger variability for low frequencies. As we previously illustrated in Fig. 4, this occurs because at the lowfrequency limit, pore fluid pressure opposition to compression is minimal, and therefore the compressibility contrast between the background and the open regions of the fractures reaches the maximal value. As a consequence, the effects of the interaction between contact areas are more important and then a change in contact areas geometry, as the one taking place in each realisation, produces a slightly different fracture compliance. As frequency increases, there is less time for fluid pressure equilibration between fracture and background, the stiffening effect of the fluid saturating the fractures increases, and the effects of the interaction between contact areas are less significant, resulting in a smaller variability of the Pwave modulus. For the three analysed frequencies, the variability of the standard deviation became approximately stable after 15 realisations, and it is lower than 1 GPa. In other words, it is always lower than 2 % of the realvalued Pwave modulus.
3.3 Comparison between realistic and simplified equivalent fracture models
In studies on the effective seismic response of fractured media, a fracture is frequently represented as a thin compliant layer of constant thickness (e.g. Schoenberg, 1980; Gurevich, 2003; Brajanovski et al., 2006; Nakagawa and Schoenberg, 2007; Barbosa et al., 2016). To analyse the validity of such simple approaches, we estimate equivalent elastic properties of a fracture from the excess compliance computed for realistic fracture models (Fig. 5). We then compare the seismic response of simple fracture models, using the derived equivalent properties, with the numerical results for realistic fractures presented in Fig. 6.
For estimating the equivalent properties of the fractures, we follow the linear slip theory (Schoenberg and Douma, 1988), where the compliance of a dry rock can be split into the contributions from the background and from the presence of fractures, and compliance is defined as the inverse of the stiffness. In this case,
where S is the total compliance matrix with coefficients S_{ij} and S_{b} is the compliance matrix of the dry background. The matrix ΔS is the excess compliance of the rock due to the presence of fractures which, in turn, is approximated as
Using Eqs. (10) and (11), and the fact that the effective stiffness matrix is the inverse of the compliance matrix, that is $\mathbf{C}={\mathbf{S}}^{\mathrm{1}}$, the normal and shear excess compliances (Z_{N} and Z_{T}), respectively, can be obtained as
where L_{m} and μ_{m} are the Pwave and shear moduli of the dry background based on properties given in Table 1. ${C}_{kk}^{\mathrm{dry}}$ values ($k=\mathrm{3},\mathrm{4}$) are numerically computed from applying corresponding compressional and shear oscillatory tests to the dry fractured models. Finally, we calculate the equivalent elastic properties from the excess compliances as (Brajanovski et al., 2005)
where ${\mathit{\mu}}_{\mathrm{fr}}^{\mathrm{eqv}}$ and ${K}_{\mathrm{fr}}^{\mathrm{eqv}}$ represent the dry shear and bulk moduli, respectively, of an equivalent fracture of constant thickness, and f_{c} is the fracture volume fraction of the models, given by the ratio between the fracture volume and the volume of the numerical volume. In this case, as the area in the plane xy of the numerical model is equal to that of the fracture (Fig. 1), f_{c} is given by the ratio between the fracture mean aperture and the height of the model (i.e. ${f}_{\mathrm{c}}={h}_{\mathrm{mean}}/H$). In order to keep the pore volume of the fractures the same for the equivalent planar fractures, we also calculate the weighted average of the fractures porosities (i.e. fracture equivalent porosity ${\mathit{\varphi}}_{\mathrm{fr}}^{\mathrm{eqv}}$) for the models shown in Fig. 5 accounting for the effect of contact areas. That is,
where ϕ_{b} is the background porosity (which is the same as for the contact areas), ϕ_{fr} is the porosity of the open regions of the fracture, ρ_{ca} is the contact area density, with ρ_{ca}=0.05 for models A and B, and ρ_{ca}=0.2 for models C and D (Fig. 5). Also, ϕ_{b}=0.1 and ϕ_{fr}=0.9 as presented in Table 1. Values obtained from Eqs. (12)–(16) are shown in Table 2.
To support the consideration of our models as realistic, the fracture compliances presented in Table 2 are in reasonable agreement with estimated values observed in field and laboratory experiments for fractured rocks (Worthington and Lubbe, 2007; Hobday and Worthington, 2012; Barbosa et al., 2019). According to the linear trend proposed by Worthington and Lubbe (2007), our models are representative of fractures with lengths ranging from a few metres to tens of metres. Furthermore, fixing a similar value for both the normal and shear compliances is usually assumed (Liu et al., 2000; Sayers, 2002); that is ${Z}_{N}/{Z}_{T}\approx \mathrm{1}$. Lubbe et al. (2008) used rock samples with artificial fractures to calculate normal and shear compliances from laboratory measurements and they found fracture compliance ratios approaching 0.5. In agreement with their results, the fractures presented in this work exhibit compliance ratios between 0.5 and 1, but closer to 0.5 (Table 2).
3.3.1 Pwave modulus analysis
After computing the equivalent dry bulk and shear moduli and porosity for each fracture (Eqs. 12–16), as well as the mean aperture, we employ the analytical solution of White et al. (1975) to quantify the Pwave modulus normal to a periodic distribution of constantthickness fractures having these equivalent fracture properties (${\mathit{\mu}}_{\mathrm{fr}}^{\mathrm{eqv}}$, ${K}_{\mathrm{fr}}^{\mathrm{eqv}}$, ${\mathit{\varphi}}_{\mathrm{fr}}^{\mathrm{eqv}}$ and h_{mean}). Figure 9 shows excellent agreement of the real part of the Pwave modulus and attenuation between fractures with variable aperture distributions (Fig. 5) and equivalent constantthickness fractures. These results show that a very simple fracture geometry, such as a thin layer of constant thickness, can approximate much more complicated fracture geometries if appropriate equivalent properties (accounting for contact area distributions) are used.
3.3.2 Stiffness matrix generalisation
To further verify and generalise the validity of the equivalent fracture model of constant thickness, we extended the methodology presented by Rubino et al. (2016) to numerically compute the effective stiffness matrix from 3D simulations for the realistic fracture models of Fig. 5. To do so, we numerically performed oscillatory relaxation tests following the methodology described in Sect. 2 but considering three normal relaxation tests (one for each direction: x, y and z) and three shear relaxation tests (shearing in the xy, xz and yz planes). We compare the results with those of the analytical solution for transverse isotropy (TI) media proposed by Krzikalla and Müller (2011). This analytical solution is based on the relaxed and unrelaxed poroelastic Backus averages of a layered porous medium consisting of a periodic distribution of a stiff background and a soft thin layer. Moreover, they showed that a single relaxation function can be used to link the relaxed and unrelaxed limits of all components of the stiffness matrix. The corresponding frequency dependence is derived from the Pwave modulus predicted by White et al. (1975). We use such soft layer to approximate a fracture of constant thickness having the equivalent properties (${\mathit{\mu}}_{\mathrm{fr}}^{\mathrm{eqv}}$, ${K}_{\mathrm{fr}}^{\mathrm{eqv}}$, ${\mathit{\varphi}}_{\mathrm{fr}}^{\mathrm{eqv}}$ and h_{mean}) obtained as described above.
Figure 10 illustrates the real part of the five independent coefficients C_{ij} for a medium with TI and C_{66}, given by ${C}_{\mathrm{66}}=\mathrm{0.5}({C}_{\mathrm{11}}{C}_{\mathrm{12}})$. The circles in the plots correspond to fractures with variable aperture distribution (Fig. 5), and the solid lines correspond to the analytical solution considering fractures of constant thickness with equivalent properties (Table 2). For brevity, fracture C was omitted due to its negligible Pwave modulus (C_{33}) dispersion. The stiffness matrix coefficients C_{11} and C_{33} dominate the stress in the medium as a response to a horizontal and vertical compression, respectively (x and z directions). C_{44} and C_{66}, on the other hand, dominate the stress as a response to a vertical and a horizontal shear deformation, respectively (yz and xy directions). The coupling coefficients C_{12} and C_{13} are also shown. The effective anisotropy of all the models is vertical transverse isotropy (VTI), although models D and B, due to their high correlation lengths, present a subtle discrepancy between C_{44} and C_{55}, which is below 0.7 GPa (omitted in the figure for brevity). Furthermore, all the stiffness coefficients of fracture models A and D are similar, which extends the conclusion about the counteracting effects of the correlation length and density of the contact areas observed for the normal Pwave modulus to the overall anisotropy of the models (Sect. 3.2). Finally, the excellent agreement between the results for models with realistic and simple geometries for all the stiffness matrix coefficients generalises our results of Sect. 3.3.1 to all incidence angles and wave modes.
Obtaining information on the hydraulic and mechanical behaviours of fractures by means of their seismic responses is an ultimate goal in fracture characterisation. PyrakNolte and Morris (2000) showed that in response to an increase in normal stress applied to a fracture, new contact areas would be created, affecting the fluid flow path through the fracture. In this work, we used similar fracture models to study the relationship between fracture microstructure and their corresponding seismic response. We showed that the effective seismic attenuation and velocity dispersion due to FPD between the fracture and the background are sensitive to changes in the correlation length and the density of contact areas. This suggests that the effects of FPD on the seismic response of fractures are potentially affected by normal stress variations and thus represent an opportunity for indirectly monitoring the associated fluid flow changes.
Milani et al. (2016) showed that in the case of aligned periodically distributed fractures having constant aperture, the models are essentially unidimensional, and no boundary effects would play a role given that the stress fields are constant within the sample as a consequence of Eq. (1). However, the distribution of the fracture apertures that we considered imply a departure from the unidimensional solution, suggesting that boundary effects may affect the results. To study these possible effects, we plotted the normalised vertical stress field in the case of the binarised model B (i) for a single repeating unit cell (RUC) and (ii) for four RUCs (not shown here). Since we found both normalised stress fields to be equal, this analysis supports the fact that no boundary effects or undesired interaction between fractures are affecting the results of the numerical relaxation tests. Moreover, this represents an extension of the results shown by Milani et al. (2016) for fractured media with 2D geometries to 3D.
We also showed that the density and correlation length of contact areas control the normal and shear fracture compliances (Fig. 6), and thus firstorder analytical solutions based on these properties, such as the model proposed by Hudson and Liu (1999), could represent a good approximation to the mechanical behaviour of realistic fractures. However, as shown in Fig. 4, the effects of the mechanical interactions between contact areas must be considered for accurately obtaining the equivalent properties of the material filling the fractures. In this sense, we showed in Fig. 10 that a simple layer of constant thickness can successfully reproduce the seismic response of fractures with intricate aperture distributions provided that appropriate elastic moduli, porosity and aperture are used. In this work, such equivalent properties are calculated from the excess compliances of the realistic fractures, and therefore they account for the effects of distribution and interaction of contact areas. Furthermore, since the linear slip model can be interpreted as an approximation of the seismic response of a constantthickness fracture (Schoenberg and Douma, 1988; Brajanovski et al., 2005; Guo et al., 2017), our results suggest that this model can also satisfactorily approximate the seismic response of the realistic fracture models considered in Fig. 5.
The aim of the present contribution was to analyse the effects of variable aperture distributions of 3D fracture models on FPD between fractures and background. To do so, we numerically quantified the effective frequencydependent stiffness matrix coefficients and seismic attenuation for realistic fracture models representing REVs of media containing periodically distributed parallel fractures. Our fracture models were characterised by aperture distributions generated using a stratified percolation algorithm and accounting for different densities and correlation lengths of contact areas. We showed that for a given density of contact areas, fractures with correlated distributions of contact areas (i.e. highest correlation length) exhibit higher Pwave modulus dispersion and seismic attenuation than those with a low correlation. Furthermore, lower Pwave modulus dispersion and seismic attenuation were observed when increasing contact area density for a given correlation length. This compensatory effect allows fractures with highly different aperture distributions to produce similar seismic responses. Moreover, although the effects of distribution of contact areas on the Pwave modulus are maximal at the lowfrequency limit, these distributions also play an important role at the highfrequency limit. We also observed that, if the distribution of contact areas is fixed, fracture mean aperture (which controls fracture volume) dominates the seismic response due to FPD effects, while the variable aperture in the open regions of the fracture has a negligible influence.
Finally, we demonstrated that a simple fracture geometry such as thin layers with constant aperture and appropriate equivalent physical properties produces the same effective anisotropic seismic response of fractures with a much more intricate geometry. The equivalent elastic properties can be obtained from the excess fracture compliances, determined according to the linear slip theory, as long as the effects of contact areas are accounted for. Our results validate the use of simple models of fractures having constant thickness for numerically simulating the effects of fractures with realistic geometries which, in turn, can significantly reduce computational cost and overcome meshing limitations.
All numerical results are reproducible by solving the equations and boundary conditions described in this work and using the fracture models provided in the Supplement. Numerical results can also be shared by contacting the first author.
To generate the fracture models shown in Fig. 5, we follow the stratified percolation approach described by Nolte and PyrakNolte (1991). To do so, we define a matrix initialised with zeros in all its cells, which represents the initial aperture distribution (aperture equal to zero) of a square fracture. The first step consists of stochastically selecting a number (x) of matrix cells. Then, each of those x cells becomes the centre of a new stochastic distribution of x cells confined to a square area around them. The squares are defined by a certain number of cells of the matrix. This process is repeated n times, called tiers, by reducing the size of the squares from one tier to the next one. It means that, if a fix number of cells (x) are stochastically located inside each of the squares corresponding to the previous tier, this finally results in selecting x^{n} cells inside the matrix. In the last tier, which is number n (corresponding to the minimum square size), a value 1 is assigned to each of the x^{n} selected cells. Moreover, as there could be overlap between the areas of the squares in all the tiers, and this is a cumulative process, at each time a value 1 is added to a cell of the matrix, the aperture of the fracture at that cell is increased. After the values of the (x) cells of the last tier (n) have been added with a 1, the contact areas are given by the matrix cells remaining with the initial zeros as values. The correlation length of the contact areas is calculated after the generation process following the methodology proposed by (Blair et al., 1993). In general, a high number of tiers (n) with a low number of cells (x) creates a correlated distribution of contact areas. On the other hand, a low number of tiers (n) with a big number of cells (x) builds a noncorrelated distribution of contact areas (PyrakNolte and Morris, 2000). The workflow for building up our models consists of importing the matrix, generated using the abovedescribed algorithm, into COMSOL Multiphysics by using an .stl format. The matrix, representing the aperture distribution of the fracture, is converted into a surface and a fracture volume is generated by symmetrically duplicating the fracture surface with respect to the horizontal middle plane. As a result of this process, we obtain a fracture volume with symmetrical walls. Finally, the obtained fracture volume is vertically scaled to yield the desired mean fracture aperture and embedded in a cubic model (Fig. 1a).
The supplement related to this article is available online at: https://doi.org/10.5194/se1013212019supplement.
This work is part of the PhD project of SL, conducted under the supervision of BQ. The paper was prepared by SL with contributions from all coauthors.
The authors declare that they have no conflict of interest.
Simón Lissa and Beatriz Quintal thank Holger Steeb and Eva Caspari for insightful discussions. J. Germán Rubino acknowledges financial support from Agencia Nacional de Promoción Científica y Tecnológica (PICT 20172976) and a visit to the University of Lausanne financed by the Fondation Herbette.
This research has been supported by the Swiss National Science Foundation (grant no. 172691).
This paper was edited by Tarje NissenMeyer and reviewed by Junxin Guo and one anonymous referee.
Amalokwu, K., Chapman, M., Best, A. I., Minshull, T. A., and Li, X. Y.: Water saturation effects on Pwave anisotropy in synthetic sandstone with aligned fractures, Geophys. J. Int., 202, 1088–1095, https://doi.org/10.1093/gji/ggv192, 2015. a
Barbosa, N. D., Rubino, J. G., Caspari, E., Milani, M., and Holliger, K.: Fluid pressure diffusion effects on the seismic reflectivity of a single fracture, J. Acoust. Soc. Am., 140, 2554–2570, https://doi.org/10.1121/1.4964339, 2016. a, b
Barbosa, N. D., Caspari, E., Rubino, J. G., Greenwood, A., Baron, L., and Holliger, K.: Estimation of fracture compliance from attenuation and velocity analysis of full waveform sonic log data, submitted, J. Geophys. Res., 124, 2738–2761, https://doi.org/10.1029/2018JB016507, 2019. a
Biot, M. A.: General theory of threedimensional consolidation, J. Appl. Phys., 12, 155–164, https://doi.org/10.1063/1.1712886, 1941. a, b, c
Biot, M. A.: Mechanics of deformation and acoustic propagation in porous media, J. Appl. Phys., 33, 1482–1498, https://doi.org/10.1063/1.1728759, 1962. a, b
Blair, S. C., Berge, P. A., and Berryman, J. G.: Twopoint correlation functions to characterize microgeometry and estimate permeabilities of synthetic and natural sandstones, Lawrence Livermore National Laboratory report, UCRLLR114182, https://doi.org/10.2172/10182383, 1993. a, b
Bourbié, T., Coussy, O., and Zinszner, B.: Acoustics of porous media, Editions Technip, Paris, 1987. a
Brajanovski, M., Gurevich, B., and Schoenberg, M.: A model for Pwave attenuation and dispersion in a porous medium permeated by aligned fractures, Geophys. J. Int., 163, 372–384, https://doi.org/10.1111/j.1365246X.2005.02722.x, 2005. a, b, c
Brajanovski, M., Müller, T. M., and Gurevich, B.: Characteristic frequencies of seismic attenuation due to waveinduced fluid flow in fractured porous media, Geophys. J. Int., 166, 574–578, https://doi.org/10.1111/j.1365246X.2006.03068.x, 2006. a
Carcione, J. M., Gurevich, B., Santos J. E., and Picotti, S.: Angular and frequencydependent wave velocity and attenuation in fractured porous media, Pure Appl. Geophys., 170, 1673–1683, https://doi.org/10.1007/s0002401206368, 2013. a
Caspari, E., Milani, M., Rubino, J. G., Müller, T. M., Quintal, B., and Holliger, K.: Numerical upscaling of frequencydependent P and Swave moduli in fractured porous media, Geophys. Prospec., 64, 369–379, https://doi.org/10.1111/13652478.12393, 2016. a
Caspari, E., Novikov, M., Lisitsa, V., Barbosa, N. D., Quintal, B., Rubino, J. G., and Holliger, K.: Attenuation mechanisms in fractured fluidsaturated porous rocks:a numerical modelling study, Geophys. Prospec., 67, 935–955, https://doi.org/10.1111/13652478.12667, 2019 a, b
Chapman, M.: Frequencydependent anisotropy due to mesoscale fractures in the presence of equant porosity, Geophys. Prospec., 51, 369–379, https://doi.org/10.1046/j.13652478.2003.00384.x, 2003. a
Guo, J., Rubino, J. G., Barbosa, N. D., Glubokovskikh, S., and Gurevich, B.: Seismic dispersion and attenuation in saturated porous rocks with aligned fractures of finite thickness: Theory and numerical simulations – Part 1: Pwave perpendicular to the fracture plane, Geophysics, 83, 49–62, https://doi.org/10.1190/geo20170065.1, 2017. a, b
Gurevich, B.: Elastic properties of saturated porous rocks with aligned fractures, J. Appl. Geophys., 54, 203–218, https://doi.org/10.1016/j.jappgeo.2002.11.002, 2003. a
Gurevich, B., Brajanovski, M., Galvin, R. J., Müller, T. M., and TomsStewart, J.: Pwave dispersion and attenuation in fractured and porous reservoirs – poroelasticity approach, Geophys. Prospec., 57, 225–237, https://doi.org/10.1111/j.13652478.2009.00785.x, 2009. a
Hobday, C. and Worthington, M.: Field measurements of normal and shear fracture compliance, Geophys. Prospec., 60, 488–499, https://doi.org/10.1111/j.13652478.2011.01000.x, 2012. a
Hudson, J. A. and Liu, E.: Effective elastic properties of heavily faulted structures, Geophysics, 64, 479–485, https://doi.org/10.1190/1.1444553, 1999. a, b, c
Hudson, J. A., Liu, E., and Crampin, S.: Transmission properties of a plane fault, Geophys. J. Int., 125, 559–566, https://doi.org/10.1111/j.1365246X.1996.tb00018.x, 1996. a
Jaeger, J. C., Cook, N. G. W., and Zimmerman, R. W.: Fundamentals of rock mechanics, 4th edn., ISBN:9780632057597, Fourth edition published 2007 by Blackwell Publishing Ltd 350 Main Street, Malden, MA 021485020, USA, 2007. a
Jänicke, R., Quintal, B., and Steeb, H.: Numerical homogenization of mesoscopic loss in poroelastic media, Eur. J. Mech. ASolid, 49, 382–395, https://doi.org/10.1016/j.euromechsol.2014.08.011, 2014. a
Krzikalla, F. and Müller, T. M.: Anisotropic PSVwave dispersion and attenuation due to interlayer flow in thinly layered porous rocks, Geophysics, 76, 135–145, https://doi.org/10.1190/1.3555077, 2011. a, b
Liu, E., Hudson, J. A., and Pointer, T.: Equivalent medium representation of fractured rock J. Geophys. Res., 105, 2981–3000, https://doi.org/10.1029/1999JB900306, 2000. a
Lubbe, R., Sothcott, J., Worthington, M. H., and McCann, C.: Laboratory estimates of normal and shear fracture compliance, Geophys. Prospec., 56, 239–247, https://doi.org/10.1111/j.13652478.2007.00688.x, 2008. a
Masson, Y. J. and Pride, S. R.: Poroelastic finite difference modeling of seismic attenuation and dispersion due to mesoscopicscale heterogeneity, J. Geophys. Res., 112, B03204, https://doi.org/10.1029/2006JB004592, 2007. a, b
Masson, Y. J. and Pride, S. R.: On the correlation between material structure and seismic attenuation anisotropy in porous media, J. Geophys. Res., 119, 2848–2870, https://doi.org/10.1002/2013JB010798, 2014. a
Masson, Y. J. and Pride, S. R.: Mapping the mechanical properties of rocks using automated microindentation tests, J. Geophys. Res., 120, 7138–7155, https://doi.org/10.1002/2015JB012248, 2015. a
Maultzsch, S., Chapman, M., Liu, E., and Li, X. Y.: Modelling frequencydependent seismic anisotropy in fluidsaturated rock with aligned fractures: implication of fracture size estimation from anisotropic measurements, Geophys. Prospec., 51, 381–392, https://doi.org/10.1046/j.13652478.2003.00386.x, 2003. a
Metz, B., Davidson, O., De Coninck, H., Loos, M., and Meyer, L.: Ipcc, 2005: IPCC special report on carbon dioxide capture and storage, Prepared by Working Group III of the Intergovernmental Panel on Climate Change, Cambridge, United Kingdom and New York, NY, USA, 2005. a
Milani, M., Rubino, J. G., Müller, T. M., Quintal, B., Caspari, E., and Holliger, K.: Representative elementary volumes for evaluating effective seismic properties of heterogeneous poroelastic media, Geophysics, 81, D21–D33, https://doi.org/10.1190/geo20150173.1, 2016. a, b
Montemagno C. D. and PyrakNolte, L. J.: Fracture network versus single fractures: Measurement of fracture geometry with Xray tomography, Phys. Chem. Earth Pt. A, 24, 575–579, https://doi.org/10.1016/S14641895(99)000824, 1999. a, b
Müller, T. M., Gurevich, B., and Lebedev, M.: Seismic wave attenuation and dispersion resulting from waveinduced flow in porous rocks – A review, Geophysics, 75, 75–147, https://doi.org/10.1190/1.3463417, 2010. a
Nakagawa, S. and Schoenberg, M. A.: Poroelastic modeling of seismic boundary conditions across a fracture, J. Acoust. Soc. Am., 122, 831–847, https://doi.org/10.1121/1.2747206, 2007. a, b
Nolte, D. and PyrakNolte, L.: Stratified continuum percolation  Scaling geometry of hierarchical cascades, Phys. Rev. A, 44, 6320–6333, https://doi.org/10.1103/PhysRevA.44.6320, 1991. a, b, c, d, e
O'Connell, R. J. and Budiansky, B.: Measures of dissipation in viscoelastic media, Geophys. Res. Lett., 5, 5–8, https://doi.org/10.1029/GL005i001p00005, 1978. a
Pride, S. R. and Berryman, J. G.: Linear dynamics of doubleporosity and dualpermeability materials, Part I: Governing equations and acoustic attenuation, Phys. Rev. E, 68, 036603, https://doi.org/10.1103/PhysRevE.68.036603, 2003. a
Pride, S. R., Berryman, J. G., and Harris, J. M.: Seismic attenuation due to waveinduced flow, J. Geophys. Res., 109, B01201, https://doi.org/10.1029/2003JB002639, 2004. a
PyrakNolte, L. J. and Morris, J. P.: Single fractures under normal stress: The relation between fracture specific stiffness and fluid flow, Int. J. Rock Mech. Min., 37, 245–262, https://doi.org/10.1016/S13651609(99)001045, 2000. a, b, c, d, e
Quintal, B., Steeb, H., Frehner, M., and Schmalholz, S. M.: Quasistatic finiteelement modeling of seismic attenuation and dispersion due to waveinduced fluid flow in poroelastic media, J. Geophys. Res., 81, 335–344, https://doi.org/10.1029/2010JB007475, 2011. a, b, c
Quintal, B., Jänicke, R., Rubino, J. G., Steeb, H., and Holliger, K.: Sensitivity of Swave attenuation to the connectivity of fractures in fluidsaturated rocks, Geophysics, 79, 15–24, https://doi.org/10.1190/geo20130409.1, 2014. a, b, c
Quintal, B., Rubino, J. G., Caspari, E., and Holliger, K.: A simple hydromechanical approach for simulating squirttype flow, Geophysics, 81, 335–344, https://doi.org/10.1190/geo20150383.1, 2016. a
Rubino, J. G., Ravazzoli, C. L., and Santos, J. E.: Equivalent viscoelastic solids for heterogeneous fluidsaturated porous rocks, Geophysics, 74, N1–N13, https://doi.org/10.1190/1.3008544, 2009. a, b
Rubino, J. G., Guarracino, L., Müller, T. M., and Holliger, K.: Do seismic waves sense fracture connectivity?, Geophys. Res. Lett., 40, 692–696, https://doi.org/10.1002/grl.50127, 2013. a, b, c
Rubino, J. G., Müller, T. M., Milani, M., and Holliger, K.: Seismic attenuation and velocity dispersion in fractured rocks: The role played by fracture contact areas, Geophys. Prospec., 62, 1278–1296, https://doi.org/10.1111/13652478.12170, 2014. a, b
Rubino, J. G., Caspari, E., Müller, T. M., Milani, M., Barbosa, N. D., and Holliger, K.: Numerical upscaling in 2D heterogeneous poroelastic rocks: Anisotropic attenuation and dispersion of seismic waves, J. Geophys. Res., 121, 6698–6721, https://doi.org/10.1002/2016JB013165, 2016. a, b
Rubino, J. G., Caspari, E., Müller, T. M., and Holliger, K.: Fracture connectivity can reduce the velocity anisotropy of seismic waves, Geophys. J. Int., 210, 223–227, https://doi.org/10.1093/gji/ggx159, 2017. a
Sayers, C. M.: Fluiddependent shearwave splitting in fractured media, Geophys. Prospec., 50, 393–402, https://doi.org/10.1046/j.13652478.2002.00324.x, 2002. a
Schoenberg, M.: Elastic wave behavior across linear slip interfaces, J. Acoust. Soc. Am., 68, 1516–1521, https://doi.org/10.1121/1.385077, 1980. a
Schoenberg, M. and Douma, J.: Elastic wave propagation in media with parallel fractures and aligned cracks, Geophys. Prospec., 36, 571–590, https://doi.org/10.1111/j.13652478.1988.tb02181.x, 1988. a, b
Schoenberg, M. and Sayers, C. M.: Seismic anisotropy of fractured rock, Geophysics, 60, 204–211, https://doi.org/10.1190/1.1443748, 1995. a
Tester, J. W., Anderson, B. J., Batchelor, A. S., Blackwell, D. D., DiPippo, R., Drake, E. M., Garnish, J., Livesay, B., Moore, M. C., Nichols, K., Petty, S., Toksoz, M. N., Veatch, R. W., Baria, R., Augustine, C., Murphy, E., Negraru, P., and Richards, M.: Impact of enhanced geothermal systems on US energy supply in the twentyfirst century, Philos. T. R. Soc. A, 365, 1057–1094, https://doi.org/10.1098/rsta.2006.1964, 2007. a
Tillotson, P., Chapman, M., Sothcott, J., Best, A. I., and Li, X.: Pore fluid viscosity effects on P and Swave anisotropy in synthetic silicacemented sandstone with aligned fractures, Geophys. Prospec., 62, 1238–1252, https://doi.org/10.1111/13652478.12194, 2014. a
White, J. E., Mikhaylova, N. G., and Lyakhovitsky, F. M.: Low frequency seismic waves in fluidsaturated layered rocks, J. Acoust. Soc. Am., 11, 654–659, https://doi.org/10.1121/1.1995164, 1975. a, b, c
Worthington, M. H. and Lubbe, R.: The scaling of fracture compliance, Geol. Soc. Spec. Pub., 270, 73–82, https://doi.org/10.1144/GSL.SP.2007.270.01.05, 2007. a, b
Zhao, L., Yao, Q., Han, D., Yan, F., and Nasser, M.: Characterizing the effect of elastic interactions on the effective elastic properties of porous, cracked rocks, Geophys. Prospec., 64, 157–169, https://doi.org/10.1111/13652478.12243, 2016. a, b
Zimmerman, R. and Main, I.: Hydromechanical behavior of fractured rocks, Int. Geophys., 89, 363–432, https://doi.org/10.1016/S00746142(03)800232, 2004. a