\documentclass{elsart}
\usepackage{subeqnarray,amsmath,amsfonts,graphicx}
\usepackage{oldgerm}
\usepackage{rotating}
\bibliographystyle{elsart-num}
\oddsidemargin -.5cm
\journal{Int. J. Rock Mech. Mining Sciences}

\renewcommand\tilde{\widetilde}
\renewcommand\a{\alpha}
\def\b{\beta}
\newcommand\bA{\mathbf{A}}
\def\bmu{\boldsymbol\mu}
\renewcommand\d{\delta}
\newcommand\D{\Delta}

\def\dfrac{\displaystyle\frac}
\def\Sum{\displaystyle\sum}
\def\ai{\left\langle\langle}
\def\ad{\right\rangle\rangle}
\def\b{\beta}
\def\k{\kappa}
\def\l{\ell}
\def\D{\Delta}
\def\d{\delta}
\def\e{\varepsilon}
\def\dint{\displaystyle\int}
\def\g{\gamma}
\def\G{\Gamma}
\def\h{{\mathbf h}}
\def\hf{ f}
\def\hR{\widehat R}
\def\hu{ u}
\def\hz{\widehat z}
\def\hZ{\widehat Z}
\def\Im{\operatorname{Im}}
\def\v{{\mathbf v}}
\def\H{{\mathbf H}}
\def\n{{\mathbf n}}
\def\R{{\mathbb R}}
\def\b{\beta}
\def\lla{\langle\langle}
\def\Lam{\Lambda}
\def\lam{\lambda}
\def\some{{x,y,z}}
\def\Lip{\operatorname{Lip}}
\def\tr{\operatorname{tr}}
\def\curl{{\nabla\times}}
\def\div{{\nabla\cdot}}
\def\opcurl{\operatorname{curl}}
\def\opdiv{\operatorname{div}}
\def\nab{{\nabla}}
\def\grad{{\nabla}}
\def\g{\gamma}
\def\q{\quad}
\def\2q{\quad\quad}
\def\3q{\quad\quad\quad}
\def\s{\sigma}
\def\sG{\mathcal{G}}
\def\sL{\mathcal{L}}
\def\sB{\mathcal{B}}
\def\sA{\mathcal{A}}
\def\cV{\mathcal{V}}
\def\sm{\setminus}
\def\NC{\mathcal{NC}}
\def\sG{\mathcal{G}}
\def\sm{\setminus}
\def\t{\theta}
\def\T{\Theta}
\def\Tau{\mathcal T}
\def\talpha{\widetilde\alpha}
\def\tA{\tilde A}
\def\tB{\tilde B}
\def\tC{\tilde C}
\def\V{\Vert}
\def\z{\zeta}
\def\xo{(x,\omega)}
\def\<{\left\langle}
\def\>{\right\rangle}
\def\div{\nabla\cdot}
\def\grad{\nabla}
\def\lma{\left\langle}
\def\la{\left\langle}
\def\lla{\langle\langle}
\def\rra{\rangle\rangle} 
\def\rma{\right\rangle} 
\def\ra{\right\rangle}
\def\qq{\quad\quad}
\def\dt{{\varDelta t}}
\def\dx{{\varDelta x}}
\def\dy{{\varDelta y}}
\def\dz{{\varDelta z}}
\def\BbbR{\mathbb R}
\def\O{\Omega}
\def\o{\omega}
\def\Oh{{\mathcal O}}
\def\mD{{\mathcal D}}
\def\mS{{\mathcal S}}
\def\mT{{\mathcal T}}
\def\mC{{\mathcal C}}
\def\mN{{\mathcal N}}
\def\p{\partial}
\def\for{\,\forall}
\def\a{\alpha}
\def\nb{\nonumber}
\def\var{\varphi}
\def\w{\omega}
\def\G{\Gamma}
\def\R{{\mathbb R}}
\def\:{{\,:\,}}
\def\Set#1{{\left\{#1\right\}}}
\def\norm#1{{\left\|#1\right\|}}
\def\qnorm#1{{\left|\left|\left|#1\right|\right|\right|}}
\def\abs#1{{\left|#1\right|}}
\def\th#1{\theta_{jk}^{#1}}
\renewcommand{\hat}{\widehat}
\newcommand{\pp}[2]{\frac{\partial {#1}}{\partial {\nu_{#2}}}} 
\newcommand{\pu}[2]{P_{\tau}u_{#1}^{#2}}
\renewcommand{\Re}{\operatorname{Re}}
\renewcommand{\Im}{\operatorname{Im}}
\newcommand{\thmrefs}[2]{Theorems~\ref{#1},~\ref{#2}}
\newcommand{\thmref}[1]{Theorem~\ref{#1}}
\newcommand{\secref}[1]{\S\ref{#1}}
\newcommand{\lemref}[1]{Lemma~\ref{#1}}
\newcommand{\rmkref}[1]{Remark~\ref{#1}}
\newcommand{\figref}[1]{Figure~\ref{#1}}
\newcommand{\ii}{{\rm i}}
\newcommand{\bsig}{{\boldsymbol \sigma}}
\newcommand{\btau}{{\boldsymbol \tau}}
\newcommand{\bvarep}{{\boldsymbol \varepsilon}}
\def\bu{{\mathbf u}}
\def\be{{\mathbf e}}

\begin{document}

\begin{frontmatter}
\title{\Large
Seismic velocity and Q anisotropy in fractured poroelastic media}
\author[igpuba,lp,pu]{Juan E. Santos\corauthref{cor}}
\author[2]{Robiel Martinez Corredor}
\author[3]{Jos\'e M. Carcione}
%\author[3]{Robiel Martinez Corredor}
\address[igpuba]{CONICET, Instituto del Gas y del Petr\'oleo, 
Facultad de Ingenier\'\i  a, Universidad de Buenos Aires,
Av. Las Heras 2214 Piso 3  C1127AAR Buenos Aires, Argentina.}
\address[lp]{Universidad
Nacional de La Plata, La Plata, Argentina.}
\address[pu]{Department of Mathematics, Purdue University, 150 N. University  Street, West Lafayette, Indiana, 47907-2067, USA.}
\address[2]{Facultad de Ingenier\'\i a, Universidad Nacional de La Plata, Argentina.}
\address[3]{Istituto Nazionale di Oceanografia e di Geofisica Sperimentale (OGS), 
Borgo Grotta Gigante 42c, 34010 Sgonico, Trieste, Italy.}
%\address[3]{Facultad de Ingenier\'\i a, Universidad Nacional de La Plata, Argentina.}
\corauth[cor]{Corresponding author. Phone: 541145143026.
\newline
 E-mail address: santos@math.purdue.edu (J. E. Santos).}

\end{frontmatter}
\newpage
\begin{frontmatter}
\begin{abstract}
A fluid-saturated poroelastic isotropic 
medium with  aligned fractures behaves as a  transversely isotropic 
and viscoelastic (TIV) 
medium when the 
predominant  wavelength is much larger
than the average distance between fractures.
A planar fracture embedded in a  fluid saturated poroelastic
 background medium   can be modeled as a extremely thin and  compliant 
porous layer. 
P-waves traveling in this type of medium 
induce fluid-pressure gradients at fractures and mesoscopic-scale heterogeneities, 
generating fluid flow and slow (diffusion) Biot waves, 
causing attenuation and dispersion of the fast modes (mesoscopic loss).
A poroelastic medium with embedded aligned fractures exhibits 
significant attenuation and dispersion effects due to this mechanism, 
which can properly be represented at the macroscale with 
an equivalent TIV medium.
In this work, we apply a set of compressibility 
and shear harmonic finite-element (FE) experiments on fractured highly heterogeneous poroelastic 
samples to determine the five 
complex and frequency dependent stiffnesses characterizing 
the equivalent medium. 
The experiments consider 
brine or patchy brine-CO$_2$ saturated samples and 
a brine saturated sample with an heterogeneous (fractal) skeleton  with  fractures. 
We show that fractures induces strong seismic velocity and Q anisotropy, both for 
qP and  qSV waves, enhanced  either by patchy saturation
or frame heterogeneity.

\end{abstract}

\end{frontmatter}

{\bf Keywords}: Fractures, poroelasticity, anisotropy, velocity dispersion, attenuation, finite elements.

\section{Introduction}
Seismic wave propagation through fractures  is an important subject in hydrocarbon 
exploration geophysics, mining and reservoir  characterization and production \cite{SD_88}. 
In particular, naturally fractured reservoirs 
have received interest in recent years, since, generally, natural fractures control the permeability of the reservoir.
In geophysical prospecting and reservoir development, knowledge of 
fracture orientation, densities and sizes is   essential 
since these factors control 
hydrocarbon production \cite{GU_03,grechka_3}. This is also important in CO$_2$ 
storage in geological formations to monitor the injected plumes 
as faults and fractures are generated, 
where  CO$_2$ can leak to the surface \cite{PIC_12}.
Among papers  presenting numerical approaches to determine 
effective media corresponding to fractured rocks,  
Grechka and Kachanov \cite{grechka_1,grechka_2} 
performed  3D static FE simulations,
summing up the individual contributions of the fractures and 
ignoring their interactions.  Also, Saenger et al. 
\cite{saenger_2} presented  numerical simulations 
in 2D and  3D media
saturated with fluids to analyze Biot's predictions in 
the high and low frequency limits of poroelasticity, while  Wenzlau et al. 
\cite{wenzlau_2010} performed FE simulations to analyze anisotropic dispersion 
and attenuation  in poroelastic materials.

A dense set of horizontal  fractures in a fluid-saturated poroelastic medium    
behaves as a    TIV medium when the average fracture distance 
is much smaller than the predominant  wavelength of the traveling waves. 
This leads to frequency and angular variations of velocity and 
attenuation of seismic waves. 
A major cause of attenuation in porous media 
is wave-induced fluid flow, which can take place
at  mesoscopic-scale heterogeneities, when the fast P-wave is converted 
into diffusion-type Biot slow waves.  Wave anelasticity and anisotropy are significant
in fractured poroelastic rocks due to  this mechanism \cite{CGSP_13}.

White et al. \cite{white_75} were  
the first to introduce the  mesoscopic-loss mechanism in the framework of Biot theory considering porous and thin plane layers. Next, Gelinsky and  Shapiro \cite{GS_97} obtained the 
relaxed and unrelaxed 
stiffnesses of the equivalent  poro-viscoelastic medium to a finely layered 
horizontally homogeneous material. 
Krzikalla and  M\"uller \cite{KM_11}
combined the two previous models assuming that 
fluid flow is perpendicular to the layering plane and independent 
of the loading direction; they obtained the five complex and frequency-dependent stiffnesses of the equivalent TIV medium \cite{CSPJGI_11}.

A planar fracture embedded in a fluid-saturated poroelastic background is  a  particular case of the thin layer problem, when one 
of the layers is very thin and compliant. FE harmonic compressibility and shear tests were first   presented in \cite{SRR_09} to obtain  a viscoelastic medium long-wavelenght equivalent to a highly heterogeneous isotropic sample.  Then, in 
\cite{PCS_10},  \cite{SCP_11} and  \cite{CSPJGI_11} the procedure was extended to determine 
long-wave equivalent media to  finely layered viscoelastic and poroelastic 
materials. The procedure used here was validated  using the analytical solution presented in \cite{KM_11}  and  \cite{CGSP_13}.

In this work, we apply  these harmonic  FE harmonic tests to saturated isotropic 
poroelastic samples having a dense set 
of horizontal  fractures modeled as very thin layers. 
The samples contained  mesoscopic-scale heterogeneities due to  patchy brine-CO$_2$ 
saturation and fractal porosity and consequently, fractal permeability and frame properties. We analyze
attenuation and velocity dispersion as function of 
frequency and propagation angle. 
 
\section{The Biot model,  equivalent medium and seismic properties}

Let us consider isotropic fluid-saturated poroelastic layers and  
let ${\mathbf u^s} ({\mathbf x}) = (u^s_1, u^s_2, u^s_3)$ 
and ${\mathbf u^f} ({\mathbf x}) = (u^f_1, u^f_2, u^f_3) $ 
indicate the  time Fourier transform of the displacement vector of the solid and 
fluid relative to the solid   frame, respectively. Here,  
if   ${\mathbf U^f}$ denotes  the fluid displacement vector, 
${\mathbf u^f} = \phi ({\mathbf U^f} -{\mathbf  u^s})$, where $\phi$ 
is the porosity.

Set ${\bu} = ({\bu^s},{\bu^f})$ and let ${\bsig}({\bu})$ and $p_f({\bu})$ denote the 
time Fourier transform of the total stress and the fluid pressure, 
respectively, and let $\be(\bu^s)$ be the  strain 
tensor of the solid phase. 
On each plane layer $n$ in a 
sequence of $N$ layers, 
the frequency-domain stress-strain relations  are \cite{carcione_07}
\begin{eqnarray}\label{st_1}
&\sigma_{kl}(\bu) =  2 \mu\, e_{kl}(\bu^s)+\delta_{kl}
\left( \lambda_{_G} \,\nabla\cdot \bu^s  + \alpha M  \nabla\cdot \bu^f \right),\\
&p_f(\bu) =  - \alpha M \nabla\cdot \bu^s -  M  \nabla\cdot \bu^f.\label{st_7}
\end{eqnarray} 
The coefficient $\mu$ is the shear modulus of the bulk 
material, considered to be equal to the shear modulus of the dry matrix. 
The other  coefficients in \eqref{st_1}-\eqref{st_7} 
can be obtained from the relations \cite{carcione_07}   
\begin{eqnarray} \label{mod.3}
&&\lambda_G= K_G - \frac{2}{3} \mu,  \quad K_G = K_m + \alpha^2 M,\\
&&\alpha = 1 - \dfrac{K_m}{K_s}, \quad  M = \left(\dfrac{\alpha- \phi}{K_s} 
+ \dfrac{\phi}{K_f}\right)^{-1},\nonumber
\end{eqnarray} 
where $K_s, K_m$ and $ K_f$  denote the bulk moduli of the solid grains,
dry matrix and saturant fluid, respectively.

Denoting by  $\o=2\pi f$  the angular frequency, 
Biot's equations of motion in the diffusive range, stated in the space-frequency
domain, are 
\begin{eqnarray}\label{motion}
&&\nabla\cdot\bsig(\bu)  = 0, \\
&& \dfrac{\ii \o \eta}{\kappa} \bu^f  + \nabla p_f(\bu)=0,\label{motion2}
\end{eqnarray}
where $\eta$ is the fluid viscosity and  $\kappa$ is the frame permeability.

Let us consider $x_1$ and $x_3$ as the horizontal and vertical coordinates,
respectively.  Gelinsky and Shapiro \cite{GS_97} showed that 
the medium behaves as   
a TI medium with the vertical symmetry axis  at long wavelengths.
They obtained the relaxed and unrelaxed limits, i.e., the low- and high-frequency 
limit real-valued stiffnesses, respectively. At all frequencies, the medium behaves as 
an equivalent  TIV  medium with 
complex and frequency-dependent 
stiffnesses, $p_{IJ}$, $I,J = 1, \ldots, 6$. For the case of 
flow normal to the  fracture layering  and 
independent of the loading direction, 
these complex stiffnesses  can be determined as  presented by 
Krzikalla and  M\"uller in \cite{KM_11} and Carcione at al. \cite{CSPJGI_11}.

Denoting by $\btau$ the stress tensor of the equivalent TIV medium and by 
$\bvarep$ 
the solid strain tensor at the macroscale, the
corresponding  stress-strain relations, 
stated in the space-frequency domain, are \cite{jose_92,carcione_07}
\begin{eqnarray}                                                           
&&\tau_{11}(u) = {p}_{11} ~\varepsilon_{11}(u^s) + p_{12} ~\varepsilon_{22}(u^s) 
+ p_{13} ~\varepsilon_{33}(u^s), \label{st_GL_1}\\
&&\tau_{22}(u) = p_{12} ~\varepsilon_{11}(u^s) + p_{11} ~\varepsilon_{22}(u^s)    
 + p_{13} ~\varepsilon_{33}(u^s), \label{st_GL_2}\\  
&&\tau_{33}(u) =  p_{13} ~\varepsilon_{11}(u^s) + p_{13} ~\varepsilon_{22}(u^s)      
+ p_{33} ~\varepsilon_{33}(u^s), \label{st_GL_3}\\   
&&\tau_{23}(u) =  2~  p_{55} ~\varepsilon_{23}(u^s), \label{st_GL_4}\\                    
&&\tau_{13}(u) = 2 ~  p_{55} ~\varepsilon_{13}(u^s),\label{st_GL_5} \\                    
&&\tau_{12}(u) = 2 ~  p_{66} ~\varepsilon_{12}(u^s). \label{st_GL_6}  
\end{eqnarray}
Here, we have assumed a closed system, for which the variation of
fluid content $\zeta = - \nabla\cdot\ \bu^f $ is equal to zero. 
This formulation  provides the complex velocities 
of the fast modes at the macroscale 
and takes into account interlayer flow effects.

The  coefficients $p_{IJ}$ in \eqref{st_GL_1}-\eqref{st_GL_6} can be determined by
applying   five compressibility and shear harmonic   
 FE  tests to a representative 2D sample of the fractured poroelastic 
material.  These tests are associated with boundary value problems for  
Biot's equations \eqref{motion} stated in the space-frequency domain. The  
 different  boundary conditions 
represent  the following virtual experiments \cite{CSPJGI_11}:

\begin{enumerate}

\item 
A compressibility test  in the  parallel  direction to the fracture layering to determine $p_{11}$.

\item A  compressibility test  in the  normal directions to the fracture layering to determine $p_{33}$.

\item
A test applying simultaneous compressions in both, the normal and parallel 
directions  to the fracture layering to determine $p_{13}$.

\item
A shear test  applied in the $(x_1, x_3)$-plane to determine $p_{55}$.

\item
A shear test in the  $(x_1, x_2)$ plane to determine $p_{66}$.

\end{enumerate}

Regarding the spatial discretization, the computational domain was 
partitioned uniformly into square cells of side length $h$. 
  The FE spaces employed to represent each component of 
the solid displacement vector ${\bu}^s$  
are locally bilinear functions which are globally continuous. The local degrees 
of freedom (DOF's)  are  the values of the components of 
${\bu}^s$ at the four corners 
of the computational cells. 

On the other hand, the relative fluid displacement 
${\bu}^f$ was represented using the vector part of the Raviart-Thomas FE 
space of zero order \cite{RT_75}.
The local DOF's  are the values of the normal component 
of  ${\bu}^f$ at the mid points 
of the faces of the computational cells. 

The arguments  presented in \cite{SRR_09} 
can be applied here to show that the error of the 
FE procedure is of the  order of 
$h^{1/2}$ in the energy norm and of the order $h$ in the $L^2$-norm. 

For a detailed description of the FE tests used in this work, we refer to   
\cite{CSPJGI_11}, where the 
model for these stiffnesses proposed by \cite{KM_11} was employed 
 to validate  the procedure.   

The complex velocities  of the equivalent TIV anisotropic medium are 
\cite{carcione_07}
\begin{eqnarray*}\label{vpPS}
%\begin{array}{l}
&&v_{\rm qP} =  (2 \bar \rho )^{-1/2} \sqrt{ p_{11} l_1^2 + p_{33} l_3^2 + p_{55} + A }, \\
&&v_{\rm qSV} =  (2 \bar \rho )^{-1/2} \sqrt{ p_{11} l_1^2 + p_{33} l_3^2 + p_{55} - A }, \\
&&v_{\rm SH} =   \bar \rho^{-1/2}\sqrt{p_{66} l_1^2 + p_{55} l_3^2}, \\
&&\hskip-.6cm A = \sqrt{ [(p_{11} - p_{55}) l_1^2 + (p_{55} - p_{33}) l_3^2]^2 + 4 [(p_{13} + p_{55}) l_1 l_3]^2},
\end{eqnarray*}
%\end{equation}
where $\bar \rho = \langle \rho \rangle$ is the 
thickness weighted average of the bulk density,
$l_1 = \sin \theta$ and $l_3 = \cos \theta$ are the directions
cosines, $\theta$ is the propagation angle between 
the wavenumber vector and the $x_3$-symmetry axis 
and the three velocities correspond to the qP, qS and SH waves, respectively.
The seismic phase velocity and quality factors  are given by
\begin{equation} \label{vph}
v_p = \left[{\rm Re} \left( \frac { 1 }{ v } \right) \right]^{-1} \ \ \ {\rm and} \quad Q = \frac { {\rm Re} (v^2) }{ {\rm Im} (v^2) }, 
\end{equation}
where $v$ represents either $v_{\rm qP}$, $v_{\rm qSV}$ or   $v_{\rm SH}$.

The energy-velocity vector ${\bf v}_e$ of the qP and qSV waves is 
\begin{equation} \label{57n}
\frac{{\bf v}_e}{v_p} = ( l_1 + l_3 \cot \psi )^{-1} \hat {\bf e}_1 +
( l_1 \tan \psi + l_3 )^{-1} \hat {\bf e}_3, 
\end{equation}
with $\psi$ being the angle between the energy-velocity vector 
and the $x_3$-axis \cite{carcione_07},  
while the energy velocity of the SH wave is \cite{carcione_07}
\begin{equation}\label{evelSHav}
{\bf v}_e =
\frac{1}{\bar \rho v_p} \left( l_1 p_{66} \hat {\bf e}_1+
l_3 p_{55} \hat {\bf e}_3 \right).
\end{equation}

\section{Numerical Results}

The FE procedures described above were implemented in FORTRAN language and 
run in the  SUN workstations of  
the Department of Mathematics at  Purdue University. 
This approach yields the five complex stiffnesses 
$p_{IJ}$ as a function  of frequency  
and the corresponding  phase velocities and  dissipation coefficients. 
For each frequency, the five discrete problems associated with the harmonic 
compressibility and shear tests    were solved using a 
public domain sparse matrix solver package.
This approach yields directly the frequency dependent  
velocities and dissipation  coefficients,  instead of 
solving Biot's equation in the space-time domain and using   
Fourier transforms to obtain the desired 
frequency domain characterization at the macroscale.

In all the experiments the numerical samples were discretized using a   
$160\times 160$ 
 uniform mesh representing 10 periods  
of 15 cm background sandstone and  1 cm  fracture thickness.
Both  background and fractures have  grain density $\rho_s = 2650$  kg/m$^{3}$, 
bulk modulus $K_s = 37$ GPa and shear modulus $\mu_s = 44$ GPa. 

The dry bulk and shear modulus of the samples were determined using 
the Krief model \cite{krief_90},
\begin{eqnarray}\label{krief}
\dfrac{K_m}{K_s} = \dfrac{\mu}{\mu_s} = (1 - \phi)^{3/(1- \phi)}.
\end{eqnarray}

A set of numerical examples consider the following  cases:

\begin{verse}

\item Case 1:  A brine-saturated sample with fractures.

\item Case 2:  A brine-CO$_2$ patchy saturated sample without fractures.

\item Case 3:  A brine-CO$_2$ patchy saturated sample with fractures.

\item Case 4: A brine saturated sample with a fractal frame and fractures. 
                  
\end{verse}

Porosity is $\phi = 0.25$ in the background and $\phi = 0.5$ in the fractures. Using 
\eqref{krief} we obtained   $K_m =1.17$ GPa and $\mu = 1.4$ GPa 
for the   background and $K_m = 0.58$ GPa and $\mu = 0.68$ GPa for  the fractures.

Permeability is obtained as
\cite{jose_2000} 
\begin{equation}\label{perm}
\kappa= \dfrac{ r_g^2 ~\phi^3}{45(1- \phi)^2}
\end{equation}
where $r_g$ =  20 $\mu$m is the average radius of the grains.

Although the algorithm has already been validated in \cite{CSPJGI_11},  
we include, for completeness, a comparison of the analytical and numerical solutions for  case  1.
Figures 1 and 2 show plots of the dissipation factors and 
energy velocities of qP and qSV waves  at 300 Hz, respectively, 
where it can be observed a very good match between the theoretical and 
numerical results.

\vskip1cm
$$
< < \text{Figure 1} > >
$$
\vskip1cm 
\vskip1cm
$$
< < \text{Figure 2} > >
$$
\vskip1cm 

Figures 3 and 4 show plots of the dissipation factors  at 50 Hz and 300 Hz,
for the cases 1, 2 and 3, while  Figures 5 and 6 display the corresponding 
polar plots of the qP energy velocity. Loss is negligible in case 1 along the direction of the fracture plane.
It can be observed strong velocity and Q anisotropy, 
with higher attenuation at  300 Hz  
and patchy brine-CO$_2$ saturation. Energy losses 
are much higher for angles between  60 and 90 degrees, i.e., for  waves  
traveling in the direction incident normal to the fracture layering. 

On the other hand, velocity anisotropy caused by the 
fractures in cases 1 and 3  is enhanced for the case of patchy saturation, 
with lower velocities when patches are present. The velocity behaves isotropically in case   2.

\vskip1cm
$$
< < \text{Figure 3} > >
$$
\vskip1cm 
\vskip1cm
$$
< < \text{Figure 4} > >
$$
\vskip1cm 

\vskip1cm
$$
< < \text{Figure 5} > >
$$
\vskip1cm 
\vskip1cm
$$
< < \text{Figure 6} > >
$$
\vskip1cm 



Figures 7 and 8  show the fluid pressure distribution 
at frequencies 50 Hz and 300 Hz, respectively, for  case  3  and  
compressions normal to the 
fracture layering. It can be observed that pressure gradients take their highest values at the fractures, and at 300 Hz remain  always higher than at 50 Hz.  
This explains the higher losses for qP waves at 300 Hz as compared with the 50 Hz experiment
observed in Figures 3 and 4.

\vskip1cm
$$
< < \text{Figure 7} > >
$$
\vskip1cm 
\vskip1cm
$$
< < \text{Figure 8} > >
$$
\vskip1cm 

Denoting by  $Q_{\rm P1}, Q_{\rm P2}$ and $Q_{\rm P3}$  the qP-quality factors associated 
with cases  1,   2 and 3 ,  Figure 9 shows the approximate validity of 
the commonly used approximation for the dissipation factors for qP waves
\begin{eqnarray}\label{qpcheck}
Q_{\rm P3}^{-1} = Q_{\rm P1}^{-1} + Q_{\rm P2}^{-1}
\end{eqnarray}
relating different attenuation mechanisms  for these cases at 300 Hz 
as function of the propagation angle. This approximation was also tested in 
\cite{gurevich_97} and \cite{CSP_12}. Numerical experiments have shown that this equation cannot be used as an approximation in Case 4.

\vskip1cm 
\vskip1cm
$$
< < \text{Figure 9} > >
$$
\vskip1cm 

Figures 10 and 11 show polar plots of the dissipation factors 
of  qSV waves at 50 Hz and 300 Hz, respectively,  for the three cases. 
For both frequencies, case   2 shows isotropic attenuation,  
while for a fractured sample brine 
or patchy saturated (cases  1 and 3 ), Q anisotropy is 
strong for angles between 30 and  
60 degrees, with  about a  50 \% increase in attenuation at 300 Hz with respect to 50 Hz. 

\vskip1cm
$$
< < \text{Figure 10} > >
$$
\vskip1cm 
\vskip1cm
$$
< < \text{Figure 11} > >
$$
\vskip1cm 

Figure 12 displays polar plots of the qSV energy velocities at 50 Hz, where velocity 
anisotropy is observed to be induced by fractures (cases  1 and 3 ), 
while patchy saturation , as expected, 
does not affect the anisotropic behavior of the  qSV velocities. 
On the other hand, case    2 shows isotropic velocity, 
with higher velocity values than for the fractured cases. 
The corresponding  qSV energy velocities at 300 Hz are almost identical to those at 
50 Hz and are not shown.

Figure 13 displays polar plots of the SH-wave energy  velocity at 50 Hz, where velocity 
is observed to be induced by fractures. Cases 1 and 3  are almost indistinguishable, while the velocity for case 2 is isotropic. 
SH energy velocities at 300 Hz are identical to those 50 Hz, since there is no energy loss ($p_{55}$ and $p_{66}$ are real \cite{CGSP_13}). 

\vskip1cm
$$
< < \text{Figure 12} > >
$$
\vskip1cm 

\vskip1cm
$$
< < \text{Figure 13} > >
$$
\vskip1cm 

Figures 14 and 15 analyze the variations in attenuation of qP waves for case 3 
 due to changes in 
CO$_2$ saturation, with attenuation coefficients computed at 300 Hz.
Figure 14 shows that an increase of  CO$_2$ saturation from 10$\%$ to 50$\%$  induces 
a noticeable decrease in attenuation for  angles 
close to the normal orientation of the fractures. Figure 15 displays the corresponding 
attenuation curves  for qSV waves, where the same decrease in attenuation is observed, 
but for angles between 30 and 60 degrees.

\vskip1cm
$$
< < \text{Figure 14} > >
$$
\vskip1cm 

\vskip1cm
$$
< < \text{Figure 15} > >
$$
\vskip1cm 

On the other hand, it was observed that 
qP energy velocities decrease for increasing CO$_2$ saturation, with the greater decreases for 
angles closer to the normal layering of the fractures. 
For qSV and SH waves, energy velocities 
show almost no change between the two CO$_2$ saturations. The related plots are omitted for brevity.

Next, we analyze the behavior of waves as a function of frequency in the 
range 1 Hz - 1 kHz at 10\% CO$_2$ saturation. 
Figures 16 and 17 display dissipation factors for cases  1,   2 and 3  for waves 
parallel (`11' waves ) and normal (`33' waves) to the fracture layering, respectively, while 
Figures 18 and 19 show the corresponding velocities.

Figure 16  indicate that   `11' waves for  case  1 (brine-saturated homogeneous background with fractures)  are lossless, while the cases of patchy saturation with and without fractures suffer similar attenuation, 
though there is a change from lower to higher attenuation for the 
patchy saturated case with fractures (case 3 ) at a frequency of about 40 Hz. 

On the other hand, the  curves for `33' waves displayed in Figure 17  
show much higher attenuation than those for `11' waves 
for the three cases. The case of patchy saturation with fractures (case 3 ) is the one exhibiting  the highest attenuation values  for all frequencies. Besides, the case of
brine-saturated homogeneous background with fractures (case  1)  shows a much higher attenuation 
than the non-fractured patchy saturated case (case   2) starting  
at a frequency of about 26 Hz, while below such frequency the opposite behavior is observed, though  with quite similar attenuation. 

\vskip1cm
$$
< < \text{Figure 16} > >
$$
\vskip1cm 


\vskip1cm
$$
< < \text{Figure 17} > >
$$
\vskip1cm 

Figure  18  shows that for the case of  brine-saturated homogeneous background 
with fractures  (case  1),  `11' velocities  are essentially independent of frequency.  
In the case of  patchy saturation with fractures (case 3 ),
velocities are always smaller than in case  1, i.e., 
the presence  of CO$_2$ patches induces a velocity decay in the whole 
range of frequencies analyzed. 
On the other hand, velocities  for patchy saturation with or without fractures
show a similar  increasing behavior with frequency. 

Figure 19 shows that  `33'  waves have the higher frequency dependent behavior for the case of a  patchy saturated fractured medium (case 3 ), with lower values than in the other two cases and increasing behavior after 10 Hz. When the sample is brine saturated  with fractures (case  1), velocities are higher than in case 3 . On the other hand, the patchy saturated case without fractures exhibits  the highest velocity values  with    
moderate increase with frequency. Thus, for `33' waves the presence of fractures induces a noticeable reduction of velocities normal to the fracture plane, either for brine or patchy saturation.  
  
\vskip1cm
$$
< < \text{Figure 18} > >
$$
\vskip1cm 


\vskip1cm
$$
< < \text{Figure 19} > >
$$
\vskip1cm

The last experiment (case 4) analyzes the effects of the presence of heterogeneities  
in the skeleton or frame of the fractured sample.
The binary fractal permeability   
is obtained  with the following logarithmic  relation \cite{russo1,russo2}
\begin{eqnarray}\label{log_rel}
\text{log} ~ \kappa(x,z) = \left\langle \text{log} ~\kappa \right\rangle + f(x,z)  
\end{eqnarray}
with $f(x,z)$ being the  fractal representing the spatial 
fluctuation of the permeability  field, chosen to be   
of fractal dimension $D = $ 2.2, correlation length 2 cm 
and average permeability 0.25 Darcy  in the background and 4.44 Darcy in the fractures.
Porosity was obtained using the Kozeny-Carman relation. 
Thus, the heterogeneous 
sample was constructed as a fractal perturbation of the sample in case  1. 
Figure 20 shows the Lam\'e coefficient $\lambda_G$ of the brine saturated fractal 
sample used in case 4.

\vskip1cm
$$
< < \text{Figure 20} > >
$$
\vskip1cm 
Figures 21 and 22 compare the qP and qSV dissipation factors of 
this case at 50 Hz  with those of case  1, while 
Figures 23 and 24 compares the corresponding energy velocities. 
As in the patchy saturation case,  frame heterogeneities induce a noticeable 
increase  in Q anisotropy for qP waves for angles normal to the fracture plane
and for qSV waves for angles between 30 and 60 degrees. From Figures 
 21 and 22, we conclude that  qSV wave attenuation is more affected  than 
 qP attenuation when frame heterogeneities are present.  Also, 
 from Figures 23 and 24  we see  the expected energy velocity reduction 
 in the heterogeneous case, and that velocity anisotropy is less affected 
 by  frame heterogeneities  than Q anisotropy.
 Concerning  SH velocities, they show anisotropy and a moderate reduction 
 in velocity and negligible attenuation 
 in the heterogeneous  case.  The SH energy velocity is similar 
 to that in case  1 of Figure 11 and the plot is omitted. 
 
 \vskip1cm
$$
< < \text{Figure 21} > >
$$
\vskip1cm 


\vskip1cm
$$
< < \text{Figure 22} > >
$$
\vskip1cm

 \vskip1cm
$$
< < \text{Figure 23} > >
$$
\vskip1cm 


\vskip1cm
$$
< < \text{Figure 24} > >
$$
\vskip1cm

\section{Conclusions}

We have presented a set of numerical quasi-static harmonic experiments to determine 
the complex and frequency dependent stiffnesses of a viscoelastic transversely 
isotropic homogeneous medium equivalent to 
a  fluid-saturated poroelastic material containing a dense set of planar 
fractures. The numerical simulators are based on the finite-element solution of 
Biot's equations in the diffusive range with boundary 
conditions representing compressibility and shear tests. 
The fractures are modeled as very thin highly permeable 
poroelastic layers of small frame moduli. 

The numerical experiments consider brine, patchy brine-CO$_2$ 
saturation and fractal frame heterogeneities, with and without fractures. Attenuation, velocity dispersion 
and anisotropy are analyzed as a function of both 
propagation angle and frequency. The effects of variation in CO$_2$ saturation has also
been analyzed.  

Strong velocity and attenuation anisotropy can be observed in the qP 
and qSV wave modes, with  attenuation 
enhanced when patches of CO$_2$ are present, 
and decreasing attenuation with increasing  CO$_2$ patchy saturation.
While the higher values of the dissipation factor for qP waves are found to be 
at angles close to the normal direction to the fracture layering, 
for  qSV  waves the higher values are found at angles between 
30 and 60 degrees. Both, qP and qSV waves, show a decrease in velocity when 
fractures are present, and Q anisotropy is more pronounced than velocity anisotropy.

There is no shear loss in the abscense of fractures, 
since the attenuation mechanism  is only 
affecting  the P waves. Fractures induce anisotropy and the P-S coupling generates 
losses in the shear waves  as well.

On the other hand, the horizontally polarized S-wave (the SH wave) 
is lossless because it is a pure mode in transversely isotropic media, 
although exhibits  velocity anisotropy, as expected.

The analysis of waves traveling parallel (`11')  and normal (`33')  
to the fracture plane as a function of frequency has been  performed in the 
range 1 Hz-1 kHz.  It has been observed 
lower dissipation factors for `11' than for `33' waves, as expected, 
and the attenuation peak for `11' waves at a lower frequency 
than that for `33' waves.

Regarding the velocities, `11' waves show a stronger dependency 
on frequency as compared with 
 `33' waves, and the presence of fractures induces a noticeable reduction 
in velocities in all the frequency range analyzed for `33' waves and at 
low frequencies up to about 90 Hz for  `11' waves.

In the last experiment, based on frame heterogeneities, 
it has been observed a similar effect on velocity and Q anisotropy 
than in the patchy saturation case, with enhanced Q anisotropy and energy 
velocity reduction in all waves, being 
the qSV waves the more affected by the fractal variations of the porous frame.  

\section{Aknowledgements}
 The work of J. E. 
Santos was partially funded by the GRANT PIP 112-20080100952 from CONICET, 
Argentina. The authors wish to thank 
Robert McGraw and J. Chapman Flack of the Computing Systems at Purdue University  
for their technical support in  running the numerical experiments.
 J. M. Carcione has been partially funded by the CO2Monitor project. 
  
\newpage

\begin{thebibliography}{999}
%%%%%%%%%%   1 
\bibitem{SD_88}
Schoenberg, M., and Douma, J., Elastic wave propagation in media with parallel fractures and aligned cracks, Geophys.
Prosp., 36, 571-590, 1988.

%%%%%%%% 2 
\bibitem{GU_03}
Gurevich, B., Elastic properties of saturated porous rocks 
with aligned fractures, J. Appl. Geophys., 54, 203-218, 2003

%%%%%%%%%%%%%%%%   3 
\bibitem{grechka_3}
Grechka, V., Comparison of the non-interaction and differential schemes in
predicting the effective elastic properties of fractured media: International
J.  of Fracture, 144, 181-188, 2007.

%%%%%%%%%%%% 4 
\bibitem{PIC_12}
Picotti, S., Carcione, J. M., Gei, D., Rossi, G., and  Santos, J. E.,
Seismic modeling to monitoring {CO}$_2$ geological storage. 
The Atzbach-Schwanenstadt gas field,
J.  Geophys. Res., 117, B06103, 2012.

%%%%%%%%%%%%%%% 5
\bibitem{grechka_1}
Grechka, V., and Kachanov, M., 
Effective elasticity of rocks with closely
spaced and intersecting cracks: Geophysics, 71,  D85-D91, 2006.

%%%%%%%%%%%%%%%%   6 

\bibitem{grechka_2}
Grechka, V., and Kachanov, M., 2006, Effective elasticity of fractured rocks: A
snapshot of the work in progress: Geophysics, 71,  W45-W58, 2006.

%%%%%%%%%%%%%   7
\bibitem{saenger_2}
Saenger E. H., Ciz, R., O. S. Kr\"uger, O. S., Schmalholz, S. M., Gurevich, B., and 
Shapiro, S. A., Finite-difference modeling of wave propagation on microscale: A
snapshot of the work in progress, Geophysics, 72, SM293-SM300, 2007.

%%%%%%%%%%%%%%%%%%   8 
\bibitem{wenzlau_2010}
Wenzlau, F., Altmann, J. B., and M\"uller, T. M., Anisotropic dispersion
and attenuation due to wave-induced flow: quasi-static finite element
modeling in poroelastic solids, J. Geophys. Res., 115, B07204,
doi:10.1029/2009JB006644, 2010.


\bibitem{CGSP_13}
Carcione, J. M., Gurevich, B., and Santos, J. E., and Picotti, S., Angular and
frequency dependent wave velocity and attenuation in fractured porous media, Pure 
 Appl.  Geophys., DOI 10.1007/s00024-012-0636-8, 2013.
 


%%%%%%%%%%%%%%%%%%%   9 

\bibitem{white_75} 
White, J. E., Mikhaylova, N. G., and Lyakhovitskiy, F. M.,  
Low-frequency seismic waves in fluid saturated layered rocks, 
Physics of the Solid Earth, 11, 654-659, 1975. 

%%%%%%%%%%%%%%   10 
\bibitem{GS_97}
Gelinsky, S., and Shapiro, S. A., Poroelastic Backus-averaging for
anisotropic, layered fluid and gas saturated sediments, 
Geophysics, 62, 1867-1878, 1997.

%%%%%%%%%%%%%%%   11

\bibitem{KM_11} 
Krzikalla, F. and M\"uller, T.,  
Anisotropic P-SV-wave dispersion and attenuation due to interlayer 
flow in thinly layered porous rocks, Geophysics, 76, W 135; doi:10.1190/1.3555077, 2011.

%%%%%%%%%%%%%%%  12
\bibitem{CSPJGI_11} 
Carcione, J. M., Santos, J. E. and Picotti, S., 
Anisotropic poroelasticity and wave-induced fluid flow. Harmonic 
finite-element simulation,  Geophys. J. Internat., 186, 1245-1254, 2011.

%%%%%%%%%%%%%%%   13 


\bibitem{PCS_10} 
Picotti, S., and Carcione, J. M., Santos, J. E., and  Gei, D., Q-anisotropy in finely-layered media, 
Geophys. Res. Lett., 37, L06302, 2010.

%%%%%%%%%%%%%%%%%%   13 a  (14)

\bibitem{SRR_09}
Santos, J. E., Rubino, J. G., and  Ravazzoli, C. L., 
A numerical upscaling procedure to estimate  effective bulk and shear
moduli in heterogeneous fluid-saturated porous media,
Comput. Methods Appl. Mech. Engrg., 198, 2067-2077, 2009.

%%%%%%%%%%%%   14 
\bibitem{SCP_11}
Santos, J. E., Carcione, J. M., and Picotti, S., Viscoelastic-stiffness tensor of
anisotropic media from oscillatory numerical experiments, 
Comput. Methods Appl. Mech. Engrg., 200, 896-904, 2011.

%%%%%%%%%%%%   15 

%%%%%%%%%%%   16 
\bibitem{carcione_07}
Carcione, J.M., 2007. Wave fields in real media: wave propagation in
anisotropic, anelastic, porous and electromagnetic media, in Handbook
of Geophysical Exploration, 2nd edn, Vol. 38, 515pp., eds Helbig, K. \&
Treitel, S., Elsevier, Oxford.

%%%%%%%%%%%%  17 



 \bibitem{jose_92}
Carcione, J. M., 
Anisotropic Q and velocity dispersion of finely layered media.
Geophys. Prosp., 40, 761-783, 1992.

%%%%%%%%%%%%%%%%%%  18 

\bibitem{RT_75}
Raviart, P. A., and Thomas, J. M., 
Mixed finite element method for 2nd order elliptic
problems, Mathematical Aspects of the Finite Element Methods, Lecture Notes
of Mathematics, vol. 606, Springer, 1975.

%%%%%%%%%%%%%  19  
\bibitem{krief_90}
Krief, M., Garat, J.  Stellingwerff, J. and Ventre, J., 
A petrophysical interpretation using the velocities 
of P and S waves (full waveform sonic), The Log Analyst, 31, 355-369, 1990.

%%%%%%%%%%%%%%   20 

\bibitem{jose_2000}
Carcione, J. M., Gurevich, B., and  Cavallini, F.,
A generalized Biot-Gassmann model for the acoustic 
properties of shaley sandstones, Geophys. Prosp., 48, 539-557, 2000.
 
%%%%%%%%%%%%%%%%%  21  

\bibitem{gurevich_97} 
Gurevich, B., Zyrianov, V. B, and Lopatnikov, S. L.,
Seismic attenuation in finely layered porous rocks: Effects 
of fluid flow and scattering, Geophysics, 62, 319-324, 1997.

%%%%%%%%%%%%%%%%%%%%%%%%%   22

\bibitem{CSP_12}
Carcione, J. M., Santos, J. E., and Picotti, S., 
Fracture induced anisotropic attenuation, Rock Mech. Rock Engrg., 45, 
929-942, 2012.

%%%%%%%%%%%%%%  23

\bibitem{russo1}
Russo, D., and Bouton, M., Statistical analysis of spatial variability
in unsaturated flow parameters, 
Water Res. Res.,  28, 1911-1925, 1992.

%%%%%%%%%%%%%%%  24 

\bibitem{russo2}
Russo, D., Russo, I., and Laufer, A., 
On the spatial variability of
parameters of the unsaturated
hydraulic conductivity,  Water Res. Res., 33,  945-956, 1997.

%%%%%%%%%%%%  25 









\end{thebibliography}

\newpage 

\listoffigures

\newpage
%\begin{figure}
%\begin{center}
%\includegraphics[width=9cm,height=9cm]
%{fig1_1000_over_qP_polar_co2_50Hz_relation_15-1_10percent_sat_negro.eps}
%\vspace{1cm}
%\caption{Analytical and numerical dissipation factor of the qP and qSV 
%waves at 300 Hz for case  1.}
%\includegraphics[width=9cm,height=9cm]
%{fig2_1000_over_qP_polar_co2_300Hz_relation_15-1_10percent_sat_negro.eps}
%\vspace{1cm}
%\caption{Analytical and numerical energy velocities  of the qP and qSV 
%waves at 300 Hz for case  1.}
%\end{center}
%\end{figure}
\newpage
\begin{figure}
\begin{center}
\includegraphics[width=9.5cm,height=9.5cm]
{fig1_1000_over_q_relat_1-15_300Hz_teorico_vs_numerico_negro.eps}
\vspace{1cm}



\caption{Dissipation factor of the qP and qSV waves at 300 Hz for case  1. The solid lines
indicate the theorical values.}

\includegraphics[width=9.5cm,height=9.5cm]
{fig2_V_relat_1-15_300Hz_teorico_vs_numerico_negro.eps}
\vspace{1cm}

\caption{Polar representation of the qP and qSV energy velocity vector 
at 300 Hz for case  1. The solid lines indicate the theorical values. }
\end{center}
\end{figure}



\newpage
\begin{figure}
\begin{center}
\includegraphics[width=9.5cm,height=9.5cm]
{fig3_1000_over_qP_polar_co2_50Hz_relation_15-1_10percent_sat_negro.eps}
\vspace{1cm}

\caption{Dissipation factor of the qP waves at 50 Hz for cases  1,   2 and 3 .}

\includegraphics[width=9.5cm,height=9.5cm]
{fig4_1000_over_qP_polar_co2_300Hz_relation_15-1_10percent_sat_negro.eps}
\vspace{1cm}

\caption{Dissipation factor of the qP waves
at 300 Hz for cases  1,   2 and 3 .}
\end{center}
\end{figure}


\begin{figure}
\begin{center}
\includegraphics[width=9.5cm,height=9.5cm]
{fig5_VqP_polar_co2_50Hz_relation_15-1_10percent_sat_negro.eps}
\vspace{1cm}
\caption{Polar representation of the qP energy velocity vector 
at 50 Hz for cases  1,   2 and 3 .}

\includegraphics[width=9.5cm,height=9.5cm]
{fig6_VqP_polar_co2_300Hz_relation_15-1_10percent_sat_negro.eps}
\vspace{1cm}
\caption{Polar representation of the qP energy velocity vector 
at 300 Hz for cases  1,   2 and 3 .}
\end{center}

\end{figure}

\newpage

\begin{figure}
\begin{center}
\includegraphics[width=13.5cm]
{fig7_presion_fluido_p33_relation_1-15_10percent_sat_negro_50Hz.eps}
\vspace{0.3cm}

\caption{Fluid pressure distribution at 50 Hz for the 
compressibility test for case 3  with compression normal to the fracture plane (`33')}

\vskip1cm 

\includegraphics[width=13.5cm]
{fig8_presion_fluido_p33_relation_1-15_10percent_sat_negro_300Hz.eps}
\vspace{0.3cm}

\caption{Fluid pressure distribution at 300 Hz for the 
compressibility test for case 3  with compression normal to the fracture plane (`33')}
\end{center}

\end{figure}

\newpage

\begin{figure}
\begin{center}
\includegraphics[width=14cm]
{fig9_1000_over_Qp3_check.eps}
\vspace{0.3cm}
\caption{Test of  the approximate validity of equation \eqref{qpcheck} for cases  1,   2 and 3 at 300 Hz.}

\end{center}

\end{figure}

\newpage

\begin{figure}
\begin{center}
\includegraphics[width=9.5cm,height=9.5cm]
{fig10_1000_over_qsv_polar_co2_50Hz_relation_15-1_10percent_sat_negro.eps}
\vspace{0.3cm}
\caption{Dissipation factor of the qSV waves  
at 50 Hz for cases  1,   2 and 3 . The lossless case 3  
is represented by a point at the origin}

\vskip1cm

\includegraphics[width=9.5cm,height=9.5cm]
{fig11_1000_over_qsv_polar_co2_300Hz_relation_15-1_10percent_sat_negro.eps}
\vspace{0.3cm}

\caption{Dissipation factor of the qSV waves
at 300 Hz for cases  1,   2 and 3 . The lossless case 3  
is represented by a point at the origin}
\end{center}

\end{figure}


\newpage

\begin{figure}
\begin{center}
\includegraphics[width=9.5cm,height=9.5cm]
{fig12_VqSV_polar_co2_50Hz_relation_15-1_10percent_sat_negro.eps}


\caption{Polar representation of the qSV energy velocity vector 
at 50 Hz for cases  1,   2 and 3 .}

\includegraphics[width=9.5cm,height=9.5cm]
{fig13_VSH_polar_co2_50Hz_relation_15-1_10percent_sat_negro.eps}

\caption{Polar representation of the SH energy velocity vector 
at 50 Hz for cases  1,   2 and 3 .}
\end{center}
\end{figure}




\newpage


\begin{figure}
\begin{center}
\includegraphics[width=9.5cm,height=9.5cm]
{fig14_1000_qP_relation_1-15_cambio_percent_sat_negro_300Hz.eps}
\vspace{0.3cm}

\caption{Dissipation factor of  qP waves for case 3 
at 300 Hz for  10\% and 50\% CO$_2$ saturation.
}

\vskip1cm

\includegraphics[width=9.5cm,height=9.5cm]
{fig15_1000_qSV_relation_1-15_cambio_percent_sat_negro_300Hz.eps}

\vskip0.3cm
\caption{Dissipation factor of  qSV  waves for case 3
at 300 Hz for  10\% and 50\% CO$_2$ saturation.}
\end{center}

\end{figure}




\newpage

\begin{figure}
\begin{center}
%\includegraphics[width=12.8cm]
\includegraphics[width=13.5cm,height=8.0cm]
{fig16_1000_over_qp11_relation_1-15_10percent_sat_negro.eps}
\vskip0.3cm
\caption{Dissipation factor of waves parallel to the fracture plane (`11' waves) 
as function of frequency}
\vskip1cm
\includegraphics[width=13.5cm,height=8.0cm]
{fig17_1000_over_qp33_relation_1-15_10percent_sat_negro.eps}
\vskip0.3cm
\caption{Dissipation factor factor of waves normal  to the fracture plane (`33'waves) 
as function of frequency}
\end{center}
\end{figure}




\newpage

\begin{figure}
\begin{center}
\includegraphics[width=13.5cm,height=8.0cm]
{fig18_vp11_relation_1-15_10percent_sat_negro.eps}
\vskip.3cm
\caption{Velocity of waves parallel to the fracture plane  (`11' waves) 
as function of frequency}

\vskip1cm
\includegraphics[width=13.5cm,height=8.0cm]
{fig19_vp33_relation_1-15_10percent_sat_negro.eps}
\vskip.3cm
\caption{Velocity of waves normal  to the fracture plane (`33'waves)
as function of frequency}
\end{center}
\end{figure}

\newpage

\begin{figure}
\begin{center}

%\includegraphics[width=12.6cm,height=10cm]
\includegraphics[width=13.5cm]
{fig20_lambda_fractal_relation_1-15_160x160_negro.eps}
\vskip.3cm
\caption{Lam\'e coefficient (GPa) for the brine-saturated 
fractal porosity-permeability sample of case 4.}
\end{center}
\end{figure}



\newpage




\begin{figure}
\begin{center}


\includegraphics[width=9.2cm,height=9.0cm]
{fig21_1000_over_qP_relation_1-15_fractal_por_perm_50Hz_negro.eps}
\vspace{0.3cm}
\caption{Dissipation factor of  qP waves  
at 50 Hz for a  brine saturated homogeneous background  with fractures (case  1)  and 
a brine saturated fractal-porosity permeability background with fractures (case 4).}

\vspace{1cm}

\includegraphics[width=9.2cm,height=9.0cm]
{fig22_1000_over_qsv_relation_1-15_fractal_por_perm_50Hz_negro.eps}
\vspace{0.3cm}
\caption{Dissipation factor of  qSV  waves  
at 50 Hz for a  brine saturated homogeneous background  with fractures (case  1)  and 
a brine saturated fractal-porosity permeability background with fractures (case 4).}
\end{center}

\end{figure}


\newpage 

\begin{figure}
\begin{center}


\includegraphics[width=9.2cm,height=9.0cm]
{fig23_VqP_relation_1-15_fractal_por_perm_50Hz_negro.eps}
\vspace{0.3cm}
\caption{Polar representation of the qP energy velocity vector 
at 50 Hz for a  brine saturated homogeneous background  with fractures (case  1)  and 
a brine saturated fractal-porosity permeability background with fractures (case 4).}

\vskip1cm 

\includegraphics[width=9.2cm,height=9.0cm]
{fig24_VqSV_relation_1-15_fractal_por_perm_50Hz_negro.eps}
\vspace{0.3cm}
\caption{Polar representation of the qSV energy velocity vector 
at 50 Hz for a  brine saturated homogeneous background  with fractures (case  1)  and 
a brine saturated fractal-porosity permeability background with fractures (case 4). }

\end{center}
\end{figure}


\end{document}


