\documentclass{elsart}
\usepackage{subeqnarray,amsmath,amsfonts,graphicx}
\usepackage{oldgerm}
\usepackage{rotating}
\bibliographystyle{elsart-num}
\oddsidemargin -.5cm
\journal{Comp. Methods Appl. Mech. Engrg.}

\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}}
\def\bv{{\mathbf v}}
\newcommand{\bnu}{{\boldsymbol \nu}}
\newcommand{\bchi}{{\boldsymbol \chi}}
\newcommand{\bpsi}{{\boldsymbol \psi}}
\newcommand{\bvarphi}{{\boldsymbol \varphi}}
\newcommand{\bgamma}{{\boldsymbol \gamma}}
\newcommand{\bTheta}{{\boldsymbol \Theta}}
\newcommand{\myfrac}[2]{\displaystyle \frac{#1}{#2}}
\begin{document}

\begin{frontmatter}
\title{\Large
Finite-element harmonic experiments to model  fractured induced anisotropy in poroelastic media}
\author[pu,igpuba,lp]{Juan E. Santos\corauthref{cor}}
\author[3]{Jos\'e M. Carcione}
\address[pu]{Department of Mathematics, Purdue University, 150 N. 
University  Street, West Lafayette, Indiana, 47907-2067, USA.}
\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[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. E-mail address: jcarcione@inogs.it}
\corauth[cor]{Corresponding author at: Department of Mathematics, Purdue University, 150 N. 
University  Street, West Lafayette, Indiana, 47907-2067, USA. Phone: 765 494 1925; Fax: 765 494 0548.
\newline
 E-mail address: santos@math.purdue.edu (J. E. Santos).}

\end{frontmatter}
\newpage
\begin{frontmatter}
\begin{abstract}
Fractures in a  fluid-saturated poroelastic -- Biot -- medium can be modeled as 
very thin highly permeable and compliant layers within a porous background.
A Biot medium containing a dense set of aligned fractures behaves as an effective 
transversely isotropic 
and viscoelastic (TIV) medium at the macroscale when the 
predominant  wavelength is much larger
than the average distance between fractures.
One  important mechanism in Biot media at seismic frequencies is wave-induced fluid flow 
generated by fast compressional waves at mesoscopic-scale heterogeneities, 
generating slow diffusion-type Biot waves. 
In this work, we present and analyze a collection of time-harmonic 
finite element experiments that take into account the effects of the presence of aligned 
fractures   and interlayer fluid flow occurring at the mesoscale, 
allowing us to  determine the complex and frequency 
dependent stiffnesses of the effective TIV medium at the macroscale.  

These  numerical upscaling experiments are defined as boundary value problems 
on  representative samples of the 
fractured material, with boundary conditions associated with compressibility and shear tests, 
which are solved using the finite element (FE) method. 
The FE space chosen to discretize each component of the solid displacement vector 
is that of globally continuous 
piecewise bilinear functions, while for the fluid  phase the vector part of 
the Raviart-Thomas-Nedelec 
space of zero order is employed.
We present results on the  uniqueness of the solution of the continuous 
and discrete problems, and derive optimal 
{\it a priori}  energy error estimates.  
First, the numerical results are validated  with those of a theory valid for  
fluid flow  perpendicular to the fracture layering  and independent 
of the loading direction, so that the attenuation mechanism can 
be represented by a single relaxation function.  
Then, the methodology is applied to cases for which no  analytical solutions are available, such 
as  a fractured Biot medium  
 saturated with brine and patches of CO$_2$ and a brine saturated sample of uniform 
background and fractures with fractal variations in their petrophysical properties.

\end{abstract}

\end{frontmatter}

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

\section{Introduction}
Fractured hydrocarbon reservoirs have been the subject of interest 
in exploration and production geophysics, since   
generally, natural fractures control the permeability of the reservoir
 \cite{GU_09}. 
In geophysical prospecting and reservoir development, knowledge of 
fracture orientation, densities and sizes is   essential 
since these factors control 
hydrocarbon production \cite{GU_03}. 
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}.
One of the important mechanisms of seismic attenuation in  
fluid-saturated  porous media 
is wave-induced fluid flow, by which the fast compressional wave is 
converted to slow (diffusive) Biot waves at mesoscopic-scale 
heterogeneities, which are larger than the pore size but much smaller 
than the predominant wavelengths of the fast compressional and shear waves \cite{CP_06,carcione_07}.

White et al. \cite{white_75} were  
the first to introduce the  mesoscopic-loss mechanism in the framework of Biot 
theory considering alternating thin 
poroelastic layers along the direction perpendicular to the layering  plane \cite{CP_06}.
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 to obtain the five complex and frequency-dependent 
stiffnesses of the equivalent TIV medium. 
Their  approach assumes a 1D character of the fluid pressure equilibration 
process which generates diffusive modes from the 
fast compressional  wave, i.e., the fluid-flow direction is perpendicular to the fracture layering. 
As a consequence,  compressional waves travelling horizontally or 
vertically or shear waves will generate fluid pressure 
in such a way as to maintain its distribution. Thus,  
  the model considers only one relaxation function, corresponding to 
the symmetry-axis  compressional wave stiffness. 
These assumptions fail for heterogeneous layers, where the  
propagation of waves may depend on direction.

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, highly permeable  and compliant. 
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. 
An analysis of wave anelasticity and anisotropy  
in fractured poroelastic rocks is presented in  \cite{CGSP_13}.

FE harmonic compressibility and shear tests is 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} and  \cite{SPC_12} the procedure is extended to determine 
long-wave equivalent media to  finely layered viscoelastic materials. 
Among other works employing numerical simulations to 
analyze dispersion, attenuation and anisotropy in 
 Biot media we mention the works by Saenger et al. 
\cite{saenger_2} and Wenzlau et al. 
\cite{wenzlau_2010}.

This work  presents and analyzes  a collection of time-harmonic 
finite  element  experiments  defined on fluid-saturated isotropic 
poroelastic samples having a dense set 
of horizontal  fractures modeled as very thin layers. The experiments take into account 
both the effects fractures  and interlayer fluid flow occurring at the mesoscale.
Each experiment is defined at the continuous level as a boundary value problem defined 
in the space-frequency domain, with boundary conditions representing compressibility and shear tests
that are solved using the FE method. 
To discretize each component of the solid displacement vector we employ the space of  of globally continuous 
piecewise bilinear functions over a quasi-regular partition of the computational domain. 
For the fluid phase, the vector part of 
the Raviart-Thomas-Nedelec 
space of zero order is used \cite{rath}. First, we demonstrate the uniqueness 
of the solution of the continuous and discrete boundary value 
problems associated with each experiment and derive 
{\it a priori} error estimates in the energy norm for the FE solutions. These estimates 
are optimal for the assumed regularity of the solution.

First, we present the validation of the FE procedure by comparison 
with the analytical solutions provided in \cite{KM_11} 
for a brine saturated   homogeneous  
sample  having  a dense set of horizontal fractures. 
Then, the methodology is applied in two cases for which no 
analytical solutions  are available, namely patchy brine-CO$_2$ saturated samples and 
brine-saturated  samples with  fractures 
having fractal variations in the  petrophysical properties.
We analyze the effect of these fluid and frame heterogeneities 
on the  attenuation and dispersion   
of the quasi-compressional wave (qP wave), 
vertically polarized quasi-shear wave (qSV wave)  and horizontally 
polarized shear wave (SH wave).
 
\section{The Biot model and the   equivalent TIV medium}

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.

Let us define the differential operator  ${\mathcal L}(\bu)$ and the matrix  
${\mathcal B}\in {\mathbf R}^4$ 
as follows:
\begin{eqnarray}\label{def_Lu}
{\mathcal L}(\bu) = \left(\nabla\cdot\bsig(\bu), \nabla p_f(\bu)\right),\quad 
%\end{eqnarray}
%\begin{eqnarray}
{\mathcal B}=
\begin{pmatrix}\label{def_B}
0 I_2& 0 I_2 \\
0 I_2 &\dfrac{\eta}{\kappa} I_2
\end{pmatrix}, 
\end{eqnarray}
where $I_2$ is the $2\times 2$ identity matrix, $\eta$ is the fluid viscosity 
and  $\kappa$ is the frame permeability.
Denoting by  $\o=2\pi f$  the angular frequency, 
Biot's equations in the diffusive range, stated in the space-frequency
domain are (in the absence of external sources) 
\begin{eqnarray}\label{diffusion}
&&i \o{\mathcal B}\bu - {\mathcal L}(\bu)= 0 ,
\end{eqnarray}
where $i = \sqrt{-1}$.
Let us consider $x_1$ and $x_3$ as the horizontal and vertical coordinates,
respectively.  As shown by Gelinsky and Shapiro (1997), the medium behaves as   
a TI medium with a vertical symmetry axis (the $x_3$-axis) at long wavelengths.
They obtaine the relaxed and unrelaxed limits, i.e., the low- and high-frequency 
limit real-valued stiffnesses, respectively. 
Assuming  a 1D character of the fluid pressure equilibration  
process, Krzikalla and  M\"uller \cite{KM_11} present a model  to 
obtain the five complex and frequency-dependent 
stiffnesses  $p_{IJ}$, $I,J = 1, \ldots, 6$,  of the equivalent TIV medium 
 that is included  in  Appendix A. 

Denote by  $\tau_{ij}(\widetilde \bu_s)$ and $\epsilon_{ij}(\widetilde \bu_s)$ the stress and strain tensor components 
of the equivalent TIV medium, where $\widetilde \bu_s$ denotes the 
solid displacement vector at the macroscale. 
The corresponding  stress-strain relations, 
stated in the space-frequency domain, are \cite{carcione_92, carcione_07} 
\begin{eqnarray}                                                           
&&\tau_{11}(\widetilde \bu_s) = {p}_{11} ~\epsilon_{11}(\widetilde \bu_s) + p_{12} ~\epsilon_{22}(\widetilde \bu_s) 
+ p_{13} ~\epsilon_{33}(\widetilde \bu_s), \label{st_GL_1}\\
&&\tau_{22}(\widetilde \bu_s) = p_{12} ~\epsilon_{11}(\widetilde \bu_s) + p_{11} ~\epsilon_{22}(\widetilde \bu_s)    
 + p_{13} ~\epsilon_{33}(\widetilde \bu_s), \label{st_GL_2}\\  
&&\tau_{33}(\widetilde \bu_s) =  p_{13} ~\epsilon_{11}(\widetilde \bu_s) + p_{13} ~\epsilon_{22}(\widetilde \bu_s)      
+ p_{33} ~\epsilon_{33}(\widetilde \bu_s), \label{st_GL_3}\\   
&&\tau_{23}(\widetilde \bu_s) =  2~  p_{55} ~\epsilon_{23}(\widetilde \bu_s), \label{st_GL_4}\\                    
&&\tau_{13}(\widetilde \bu_s) = 2 ~  p_{55} ~\epsilon_{13}(\widetilde \bu_s),\label{st_GL_5} \\                    
&&\tau_{12}(\widetilde \bu_s) = 2 ~  p_{66} ~\epsilon_{12}(\widetilde \bu_s). \label{st_GL_6}  
\end{eqnarray}   
In \eqref{st_GL_1}-\eqref{st_GL_6} we have assumed a closed system.  
This can be done for the undrained composite medium, for which the variation of
fluid content $\zeta = - \nabla\cdot \bu_f $ is equal to zero. 
This approach provides the complex velocities of the fast qP, qSV and SH  modes.
The $p_{IJ}$ are the complex and frequency-dependent 
Voigt stiffnesses  
to be determined with the harmonic experiments.
In the next sections 
we present and analyze a FE procedure to determine the coefficients 
in equations (\ref{st_GL_1})-(\ref{st_GL_6}) and the
corresponding phase velocities and quality factors. 
These properties, which depend on frequency and propagation direction, are given in Appendix B.
We show that to  determine  the $p_{IJ}$'s  
it is sufficient to perform a collection of FE harmonic  tests 
on representative 2D samples of the viscoelastic material.

\section{Determination of the stiffnesses}
Set $\G = \G^L\cup\G^B\cup\G^R\cup\G^T$, where  
\begin{eqnarray}
&&\G^L = \{ (x_1, x_3)\in \G: x_1=0\}, \quad \G^R = \{ (x_1, x_3)\in \G: x_1=L\},\nonumber\\
&&\G^B = \{ (x_1, x_3)\in \G: x_3=0\},\quad \G^T = \{ (x_1, x_3)\in \G: x_3=L\}.\nonumber
\end{eqnarray}

Denote by $\bnu$ the unit outer normal on $\G$ and let   
 $\bchi$ be a   unit tangent  on $\G$  so that $\{\bnu,\bchi\}$ is 
an orthonormal system on $\G$.

i) To determine the complex  coefficient $p_{33}$  
let us consider the solution of \eqref{diffusion} in $\O$    
together with the following boundary conditions 
\begin{eqnarray}\label{bc1a}
&&\bsig(\bu) \bnu \cdot \bnu = -\Delta P, \quad (x_1, x_3)\in \G^T,\label{bc1b}\\
&&\bsig(\bu) \bnu \cdot \bchi = 0, \quad (x_1, x_3)\in \G,\label{bc1c}\\
&&\bu_s\cdot\bnu = 0, \quad (x_1, x_3)\in \G^L\cup\G^R\cup\G^B,\label{bc1e}\\
&&\bu_f\cdot\bnu=0, \quad (x_1, x_3)\in \G\label{bc1g}.
\end{eqnarray}
Denoting by $V$ the original volume of the sample, 
its (complex) oscillatory volume change, 
$\Delta V(\omega)$, allows us to define $p_{33}$ 
 by using the relation

\begin{equation}\label{K-dv}
\frac{\Delta V(\omega)}{V}=-\frac{\Delta P}{p_{33}(\omega)},
\end{equation}
valid for a viscoelastic homogeneous medium in the quasistatic case.

After solving  \eqref{diffusion} with the boundary
conditions \eqref{bc1b}-\eqref{bc1g}, the vertical displacements  $u_{s,3}(x_1,L,\omega)$ 
on  $\Gamma^T$ allow 
 us to obtain an average vertical displacement $u^{T}_{s,3}(\omega)$ suffered by the 
boundary
$\Gamma^T$. Then, for each frequency $\omega$, the volume change produced 
by the compressibility test
can be approximated by $\Delta V(\omega)\approx L u^{T}_{s,3}(\omega)$, 
which enable us to compute  $p_{33}(\omega)$ by using the relation (\ref{K-dv}). 

ii) To determine $p_{11}$,  we solve \eqref{diffusion}  
in $\O$ together with the  boundary conditions 
\begin{eqnarray}\label{bc1a_11}
&&\bsig(\bu) \bnu \cdot \bnu = -\Delta P, \quad (x_1, x_3)\in \G^R,\label{bc1b_11}\\
&&\bsig(\bu) \bnu \cdot \bchi = 0, \quad (x_1, x_3)\in \G,\label{bc1c_11}\\
%&&\bsig(\bu) \bnu \cdot \bchi = 0, \quad (x_1, x_3)\in \G^B\cup\G^T,\label{bc1d_11}\\
&&\bu_s\cdot\bnu = 0, \quad (x_1, x_3)\in \G^L\cup\G^B\cup \G^T,
\label{bc1e_11}\\
%&&\bu_s = 0, \quad (x_1, x_3)\in \G^L,\label{bc1f_11}\\
&&\bu_f\cdot\bnu=0, \quad (x_1, x_3)\in \G\label{bc1g_11}.
\end{eqnarray}
In this experiment $~\epsilon_{33}(\bu_s)= ~\epsilon_{22}(\bu_s) 
= \nabla\cdot \bu_f= 0$
and from  \eqref{st_GL_1} we see that this experiment determines $p_{11}$ 
as  indicated for $p_{33}$ measuring the oscillatory volume change.

iii) To determine $p_{55}$ 
let us consider the solution of \eqref{diffusion}   
in $\O$  with 
the following boundary conditions 
\begin{eqnarray}\label{s_bc1a}
&&-\bsig(\bu) \bnu  = {\mathbf g}, \quad (x_1, x_3)\in \G^T\cup\G^L\cup\G^R,\label{s_bc1b}\\
&&\bu_s = 0, \quad (x_1, x_3)\in \G^B,\label{s_bc1f}\\
&&\bu_f\cdot\bnu=0, \quad (x_1, x_3)\in \G\label{s_bc1g}, 
\end{eqnarray}
where 
\begin{eqnarray*}
{\mathbf g}=\begin{cases} (0, \Delta p),\quad (x_1, x_3)\in \G^L,\\
(0, -\Delta p),\quad (x_1, x_3)\in \G^R,\\
(-\Delta p,0), \quad (x_1, x_3)\in \G^T.\end{cases} 
\end{eqnarray*}
The change in shape of the rock sample allows to obtain   
$p_{55}(\omega)$ by using the relation
\begin{equation}\label{tg}
\text{tg}(\beta\omega))=\frac{\Delta p}{p_{55}(\omega)},
\end{equation}
where $\beta(\omega)$ is the departure angle between the original positions of the 
lateral boundaries
and those after applying the shear stresses (see, for example, \cite{kolsky63}). 
Equation \eqref{tg} holds for this experiment in a viscoelastic homogeneous media 
in the  quasistatic approximation.

The horizontal displacements 
$u_{s,1}(x_1, L,\omega)$ at the top 
boundary $\Gamma^T$ allow  us to obtain, for each frequency,  an average horizontal 
displacement  $u^T_{s,1}(\omega)$ suffered by the 
boundary $\Gamma^T$. This average value allows us to approximate the 
change in shape suffered by
the sample, given by $\text{tg}(\beta(\omega))\approx u^T_{s,1}(\omega)/L$, 
which from \eqref{tg}
let us estimate $p_{55}(\omega)$.

v) To determine $p_{13}$ we solve \eqref{diffusion} in $\O$ 
with  the boundary conditions 
\begin{eqnarray}\label{bc1a_13}
&&\bsig(\bu) \bnu \cdot \bnu = -\Delta P, \quad (x_1, x_3)\in \G^R\cup \G^T,
\label{bc1b_13}\\
&&\bsig(\bu) \bnu \cdot \bchi = 0, \quad (x_1, x_3)\in \G,\label{bc1c_13}\\
&&\bu_s\cdot\bnu = 0, \quad (x_1, x_3)\in \G^L\cup\G^B,\label{bc1e_13}\\
&&\bu_f\cdot\bnu=0, \quad (x_1, x_3)\in \G\label{bc1g_13}.
\end{eqnarray}

Thus, in this experiment $\epsilon_{22}(\bu_s)= \nabla\cdot \bu_f=0$, 
and from \eqref{st_GL_1} and 
\eqref{st_GL_3} we get 
\begin{eqnarray}\label{st_backus_2d}
&&\tau_{11} =  c_{11} \epsilon_{11} +  c_{13} \epsilon_{33} \\
&&\tau_{33} = c_{13} \epsilon_{11} +  c_{33} \epsilon_{33},\nonumber
\end{eqnarray}
where $\epsilon_{11}$ and $\epsilon_{33}$ are the (macroscale) 
strain components at the  right lateral side and  top side 
 of the sample, respectively. Then from \eqref{st_backus_2d} and the fact that 
 $\tau_{11}=\tau_{33}= - \Delta P$ (c.f.\eqref{bc1b_13}) we obtain
$ p_{13}(\o)$ as
\begin{equation} \label{p13}
p_{13}(\o) = \frac {c_{11} \epsilon_{11} - c_{33} \epsilon_{33}}{\epsilon_{11} - \epsilon_{33}}.
\end{equation} 

iv) To determine $ p_{66}$, since this stiffness is associated with shear waves traveling in 
the $(x_1, x_2)$-plane, 
we consider an homogeneous horizontal slab in the $x_2$-direction and an homogeneous sample $\O_2=(0,L)^2$ 
in the $(x_1, x_2)$-plane, 
with boundary $\G_2 = \G_2^L\cup\G_2^B\cup\G_2^R\cup\G_2^T$, 
where  
\begin{eqnarray}
&&\G_2^L = \{ (x_1, x_2)\in \G: x_1=0\}, \quad \G_2^R = \{ (x_1, x_2)\in \G: x_1=L\},
\nonumber\\
&&\G_2^B = \{ (x_1, x_2)\in \G: x_2=0\},\quad \G^T = \{ (x_1, x_2)\in \G: x_2=L\}.
\nonumber
\end{eqnarray}
Then let us consider the solution of \eqref{diffusion}  
 in $\O_2$  with  the following boundary conditions 
\begin{eqnarray}\label{s_bc1a_2}
&&-\bsig(\bu) \bnu  = {\mathbf g_2}, \quad (x_1, x_2) \in \G^T_2\cup\G^L_2\cup\G^R_2,\label{s_bc1b_2}\\
&&\bu_s = 0, \quad (x_1, x_2)\in \G^B_2,\label{s_bc1f_2}\\
&&\bu_f\cdot\bnu=0, \quad (x_1, x_2)\in \G\label{s_bc1g_2}, 
\end{eqnarray}
where 
\begin{eqnarray*}
{\mathbf g_2} = \begin{cases} (0, \Delta p),\quad (x_1, x_2)\in \G^L_2,\\
(0, -\Delta p),\quad (x_1, x_2)\in \G^R_2,\\
(-\Delta p,0), \quad (x_1, x_2)\in \G^T_2.\end{cases} 
\end{eqnarray*}
Then, we proceed as indicated for  $p_{55}(\omega)$.

The  stiffnesses coefficients $p_{IJ}$ allow us to to calculate the 
wave velocities and quality factors of the effective TIV medium as 
explained in  Appendix B.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{A variational formulation}
In order to state a variational formulation  
we need to introduce some notation. For  $X\subset \mathbb R^d$ with boundary $\partial X$, let $(\cdot,\cdot)_X$ and 
$\left\langle\cdot,\cdot\right\rangle_{\partial X}$
denote the complex $L^2(X)$ and $L^2(\partial X)$ inner products for scalar, 
vector, or matrix valued functions.
Also, for $s\in\mathbb R$,
$\|\cdot\|_{s,X}$ and $|\cdot|_{s,X}$ will denote  
the usual norm and seminorm for the Sobolev space $H^s(X)$,  \cite{adams_75}.
In addition, if $X=\O$ or $X=\G$, the subscript $X$ may be omitted such that
$(\cdot,\cdot)=(\cdot,\cdot)_\O$ or  
$\<\cdot,\cdot\>=\<\cdot,\cdot\>_\G$.

Let us introduce the following closed subspaces 
of $[H^1(\O)]^2$ and $[H^1(\O_2)]^2$: 
\begin{eqnarray*}
&&{\mathcal W}_{11}(\O)= \{\bv\in [H^1(\O)]^2: \, 
\bv\cdot\bnu = 0 \,\,\,  \text{on} \, \,\G^B\cup \G^T\cup\G^L\},
\\ \\
&&{\mathcal W}_{33}(\O)= \{\bv\in [H^1(\O)]^2: 
\bv\cdot\bnu = 0 \,\,\,  \text{on} 
\, \,\G^L\cup \G^R\cup \G^B \},
\\ \\
&&{\mathcal W}_{13}(\O)= \{\bv\in [H^1(\O)]^2: \,\bv\cdot\bnu = 0 \,\,\,  \text{on} \,
\,\G^L\cup \G^B\},
\\ \\
&&{\mathcal W}_{55}(\O)= \{\bv\in [H^1(\O)]^2: 
\,\bv=0  \,\,\,\text{on} \,\,\G^B\},
\\ \\ 
&&{\mathcal W}_{66}(\O_2)= \{\bv\in [H^1(\O_2)]^2:\,\, \bv=0  \,\,\,\text{on} \,\,\G_2^B\}.
\end{eqnarray*}
Also, let 
\begin{eqnarray*}
&& H_0(\opdiv;\O)=\{\bv\in H(\opdiv;\O):\, \bv\cdot\bnu = 0
\,\,\,\text{on} \,\,\G\},
\\ \\
&& H_0(\opdiv;\O_2)=\{\bv\in H(\opdiv;\O_2):\,  \bv\cdot\bnu = 0
\,\,\,\text{on} \,\,\G_2\},
\\ \\
&& H^1(\opdiv;\O)=\{\bv\in [H^1(\O)]^2:\, \nabla\cdot\bv
\in H^1(\O)\},
\end{eqnarray*}
and  for $(I,J) = (1, 1), (3, 3), (1, 3), (5, 5)$ let
\begin{eqnarray*}\label{cont_spaces}
&&{\mathcal Z}_{IJ}(\O) = {\mathcal W}_{IJ}(\O)\times  H_0(\opdiv;\O).
\end{eqnarray*}
Also, let 
\begin{eqnarray*}\label{cont_spaces_66}
&&{\mathcal Z}_{66}(\O_2) = {\mathcal W}_{66}(\O_2)\times  H_0(\opdiv;\O_2).
\end{eqnarray*}

To obtain  our variational formulation associated with $p_{33}$, 
multiply equation (\ref{diffusion}) by $\bv=\left(\bv_{s}, \bv_{f}\right)\in
 {\mathcal Z}_{33}(\O)$, 
use integration by parts  and apply the boundary conditions \eqref{bc1b}-\eqref{bc1g}
to see get  {\it the weak form}: 
find $\bu^{(33)} = (\bu_s^{(33)}, \bu_f^{(33)}) \in {\mathcal Z}_{33}(\O)$
such that  
\begin{eqnarray}\label{weak_c33}
&&\Lambda(\bu^{(33)}, \bv) \equiv  i \o \left(\dfrac{\eta}{\kappa} \bu_f^{(33)}, \bv_f\right)\\ 
&&+ \sum_{s,t}\left(\sigma_{st}(\bu^{(33)}), 
e_{st}(\bv_{s})\right)
 - \left(p_f(\bu^{(33)}),\nabla\cdot \bv_{f})\right)\nonumber\\
&&= -\left\langle \Delta P, \bv_{s}\cdot\bnu\right\rangle_{\G^T}, \qquad  
\forall \,\,\bv=\left(\bv_{s}, \bv_{f}\right)\in 
{\mathcal Z}_{33}(\O)\nonumber .
\end{eqnarray}
Note  that  in \eqref{weak_c33},  we can write 
\begin{eqnarray}
&&\sum_{l=1}^2\sum_{s,t}\left(\sigma_{st}(\bu^{(33)}), e_{st}(\bv^{s})
\right)
 - \left(p_f(\bu^{(33)}),\nabla\cdot \bv^{f})\right)_{R^{(l)}}\label{strain}
\\
&&\qquad  
= \left({\mathbf E} ~ \tilde \be(\bu^{(33)}), \tilde e(\bv)\right).\nonumber
\end{eqnarray}
In \eqref{strain}, the   
matrix ${\mathbf E}$ = and the column vector  $\widetilde \be((\bu^{(33)}))$
are defined  by
\begin{eqnarray}\label{def_E}
{\mathbf E}=\begin{pmatrix}
\lambda_G+2\mu&\lambda_G&\alpha M&0 \\
\lambda_G&\lambda_G+2\mu&\alpha M&0\\
\alpha M&\alpha M&M&0\\
0&0&0&4 \mu
\end{pmatrix}&,&\,\,
\widetilde \be(\bu^{(33)}) = \begin{pmatrix}
e_{11}(\bu_s^{(33)}) \\ e_{33}(\bu_s^{(33)}) \\ \nabla\cdot \bu_f^{(33)}\\ e_{13}(\bu_s^{(33)})
\end{pmatrix}.
\end{eqnarray}
Note  that  ${\mathbf E}$ is 
positive definite since it is associated with the strain energy density.

Thus we can state \eqref{weak_c33} in the equivalent form:
find $\bu^{(33)} = (\bu_s^{(33)}, \bu_f^{(33)}) \in {\mathcal Z}_{33}(\O)$
such that  
\begin{eqnarray}\label{weak_c33a}
&&\Lambda(\bu^{(33)}, \bv) =  i \o \left(\dfrac{\eta}{\kappa} \bu_f^{(33)}, \bv_f\right) +
\left({\mathbf E} ~ \tilde \be(\bu^{(33)}), \tilde \be(\bv)\right).
\end{eqnarray}

Similarly, we obtain the variational formulations for the other $p_{IJ}$'s:
\begin{itemize}
\item
$p_{(11)}$:
\newline
\noindent
Find $\bu^{(11)} = (\bu_s^{(11)}, \bu_f^{(11)}) \in {\mathcal Z}_{11}(\O)$
such that  
\begin{eqnarray}\label{weak_c11}
&&\Lambda(\bu^{(11)}, \bv) = -\left\langle \Delta P, \bv_{s}\cdot\bnu\right\rangle_{\G^R}, \qquad  
\forall \,\,\bv=\left(\bv_{s}, \bv_{f}\right)\,\,\in 
{\mathcal Z}_{11}(\O).
\end{eqnarray}

\item
$p_{(13)}$:
\newline
\noindent
Find $\bu^{(13)} = (\bu_s^{(13)}, \bu_f^{(13)}) \in {\mathcal Z}_{13}(\O)$
such that  
\begin{eqnarray}\label{weak_c13}
&&\Lambda(\bu^{(13)}, \bv) = -\left\langle \Delta P, \bv_{s}\cdot\bnu\right
\rangle_{\G^R\cup\G^T}, \qquad  
\forall \,\,\bv=\left(\bv_{s}, \bv_{f}\right)\in 
{\mathcal Z}_{13}(\O).
\end{eqnarray}

\item
$p_{(55)}$:
\newline
\noindent
Find $\bu^{(55)} = (\bu_s^{(55)}, \bu_f^{(55)}) \in {\mathcal Z}_{55}(\O)$
such that  
\begin{eqnarray}\label{weak_c55}
&&\Lambda(\bu^{(55)}, \bv) = -\left\langle {\mathbf g}, \bv_{s}\right
\rangle_{\G^R\cup\G^T}, \qquad  
\forall \,\,\bv=\left(\bv_{s}, \bv_{f}\right)\in 
{\mathcal Z}_{55}(\O).
\end{eqnarray}

\item
$p_{(66)}$:
\newline
\noindent
Find $\bu^{(66)} = (\bu_s^{(66)}, \bu_f^{(66)}) \in {\mathcal Z}_{66}(\O_2)$
such that  
\begin{eqnarray}\label{weak_c66}
&&\Lambda(\bu^{(66)}, \bv) = -\left\langle {\mathbf g_2}, \bv_{s}\right
\rangle_{\G^R\cup\G^T}, \qquad  
\forall \,\,\bv=\left(\bv_{s}, \bv_{f}\right)\in 
{\mathcal Z}_{66}(\O_2).
\end{eqnarray} 
\end{itemize}

\subsection{Uniqueness of the solution of the variational problems}

The five boundary-value 
problems   formulated above are associated with second-order operators  having  
boundary data in $L^2(\O)$. Existence will be assumed for the solution of these problems.
Also, following \cite{kenig_88} we will assume that 
 $\bu_s^{(IJ)} \in [H^{3/2}]^2$. On the other hand, it will be assumed that
 $\bu_f^{(IJ)} \in H^1(\opdiv;\O)$.  These assumptions will be used in the proof
 of the theorem stated below.
\begin{thm}\label{thm;uniqness}
Assume that  
$\bu_s^{(IJ)} \in [H^{3/2}(\O)]^2, \bu_f^{(IJ)} \in H^1(\opdiv;\O)$ 
 for $(I,J) = (1,1), (3,3), (1,3), (5, 5)$ 
 and that $\bu_s^{(66)} \in [H^{3/2}(\O_2)]^2, \bu_f^{(66)} \in H^1(\opdiv;\O_2)$ .
 Also assume  that 
the matrix $\mathbf E$ in \eqref{def_E} is positive definite. 
Then the solution of  problems \eqref{weak_c33a}, \eqref{weak_c11}, \eqref{weak_c13}, 
\eqref{weak_c55}  and  \eqref{weak_c66} is unique.
\end{thm}
{\it Proof:} 
To analyze the uniqueness of the solution of \eqref{weak_c33a}, set 
$\Delta P=0$ and  choose $\bv=\bu^{(33)}$ in \eqref{weak_c33a} to obtain the equation
\begin{eqnarray}\label{uniq_1}
i \o \left(\dfrac{\eta}{\kappa} \bu_f^{(33)}, \bu_f^{(33)}\right) +
\left({\mathbf E} ~ \tilde \be(\bu^{(33)}), \tilde \be(\bu^{(33)})\right)= 0.
\end{eqnarray}
Choose the imaginary part in \eqref{uniq_1} to conclude that
\begin{eqnarray}\label{uniq_2}
\|\bu_f^{(33)}\|_0= 0.
\end{eqnarray}
Using \eqref{uniq_2} in \eqref{uniq_1} we obtain

\begin{eqnarray}\label{uniq_1a}
\left({\mathbf {\widehat E}} ~ {\widehat  \be}(\bu_s^{(33)}), {\widehat  \be}(\bu_s^{(33)})\right)= 0,
\end{eqnarray}
where
\begin{eqnarray}\label{def_hatE}
{\mathbf {\widehat E}}=\begin{pmatrix}
\lambda_G+2\mu&\lambda_G&0 \\
\lambda_G&\lambda_G+2\mu&0\\
0&0&4\mu
\end{pmatrix}&,&\,\, \ \ \
{\widehat \be}(\bu_s^{(33)}) = \begin{pmatrix}
e_{11}(\bu_s^{(33)}) \\ e_{33}(\bu_s^{(33)}) \\ e_{13}(\bu_s^{(33)})
\end{pmatrix}.
\end{eqnarray}
Next, since the matrix $\mathbf E$ is positive, so it is  the 
matrix ${\mathbf {\widehat E}}$,  which from \eqref{uniq_1a}  allow us to conclude 
that
\begin{eqnarray}\label{uniq_4}
\|e_{11}(\bu_s^{(33)})\|_0 = \|e_{33}(\bu_s^{(33)})\|_0 = \|e_{13}(\bu_s^{(33)})\|_0 = 0.
\end{eqnarray}
Now  use  the Sobolev embedding  \cite{adams_75} 
\begin{eqnarray}\label{sobolev}
H^{3/2}(\O)\rightarrow C^0(\overline \O)
\end{eqnarray}
and \eqref{uniq_4}  to see that for some constants $A, B, C$ 
\begin{eqnarray}\label{uniq_5}
u_{s,1}^{(33)}(x_1, x_3) = C x_3 + B, \, \,  u_{s,3}^{(33)}(x_1, x_3) = -C x_1 + B, 
\, \forall (x_1, x_3)\in \overline \O. 
\end{eqnarray}
Next using the boundary condition \eqref{bc1e} we see that $A= B = C = 0$, so that 
\begin{eqnarray}\label{uniq_6}
\bu_s^{(33)} = 0, \quad \in \overline\O.
\end{eqnarray}
Combining \eqref{uniq_2} and \eqref{uniq_5}   we conclude that 
uniqueness holds for the solution of \eqref{weak_c33a}. 
Uniqueness for the solution of \eqref{weak_c11} and \eqref{weak_c13} 
follows with the same argument.

Let us analyze  uniqueness for  the solution of \eqref{weak_c55}. Set  
$\Delta P=0$ and  choose $\bv=\bu^{(55)}$ in \eqref{weak_c55}.  Repeating the argument 
given to show uniqueness  for $\bv=\bu^{(33)}$ we obtain 
\begin{eqnarray}\label{uniq_7}
&&\|\bu_f^{(55)}\|_0 = 0,\\
&&\|e_{11}(\bu_s^{(55)})\|_0 = \|e_{33}(\bu_s^{(55)})\|_0 = \|e_{13}(\bu_s^{(55)})\|_0 = 0.\nonumber
\end{eqnarray}
Next, note that since elements in ${\mathcal Z}_{55}(\O)$  vanish on $\G^B$,  
\begin{eqnarray}\label{equiv}
\|| \bv \|| = \left( \sum_{k,l} \int_{\O} |\e_{kl}(\bv)|^2 d\O  \right)^{1/2} 
\end{eqnarray}
defines a norm on ${\mathcal Z}_{55}(\O)$ equivalent to the $H^1$-norm in $\O$, \cite{ciarlet80}. Thus from \eqref{uniq_7} we see that
\begin{eqnarray}\label{uniq_8}
\|\bu_s^{(55)}\|_1 = 0.
\end{eqnarray}
>From \eqref{uniq_7} and \eqref{uniq_8} we see that uniqueness holds for the solution 
of \eqref{weak_c55}.
Uniqueness for the solution of \eqref{weak_c66} follows with identical argument. 
This completes the proof.

\section{The finite element method}
Let ${\mathcal T}^h(\O)$ be a  non-overlapping  partition of $\O$ 
into rectangles $\O_j$ of diameter bounded
by $h$ such that $\overline \O = \cup_J^J \overline\O_j$.
Denote by  $\G_{jk} = \p\O_j\cap\p\O_k$  the common side of two adjacent 
rectangles $\O_j$ and $\O_k$. Also, let $\G_j = \p\O_j\cap \G$.

We employ the space of globally continuous piecewise bilinear polynomials, 
denoted below below by  
${\mathcal W}^{h}_{33}(\O)$ to approximate each component of the 
solid displacement $\bu_s$, while the vector part of 
the  Raviart-Thomas-Nedelec space  
${\mathcal V}^h_{33}(\O)$ of zero order is used to approximate 
the fluid displacement vector $\bu_f$ \cite{rath}. More specifically, let 
\begin{eqnarray*}
&&{\mathcal W}^{h}_{11}(\O)= \{ \bv_s: \bv_s|_{\O_j}\in [P_{1,1}(\O_j)]^2, \,\, 
\bv_s\cdot\bnu = 0 \,\,\text{on}\,\, \G^B\cup\G^T\G^L\}\cap [C^0(\overline \O)]^2,\\ \\
&&{\mathcal W}^{h}_{33}(\O)= \{ \bv_s: \bv_s|_{\O_j}\in [P_{1,1}(\O_j)]^2, \,\, 
\bv_s\cdot\bnu = 0  \,\,\text{on}\,\,\G^L\cup\G^R\G^B\}\cap [C^0(\overline \O)]^2,\\ \\
&&{\mathcal W}^{h}_{13}(\O)= \{ \bv_s: \bv_s|_{\O_j}\in [P_{1,1}(\O_j)]^2, \,\, 
\bv_s\cdot\bnu = 0  \,\,\text{on}\,\,\G^L\cup\G^B\}\cap [C^0(\overline \O)]^2,\\ \\
&&{\mathcal W}^{h}_{55}(\O)= \{ \bv_s: \bv_s|_{\O_j}\in [P_{1,1}(\O_j)]^2, \,\, 
\bv_s\cdot\bnu = 0  \,\,\text{on}\,\, \G^B \}\cap [C^0(\overline \O)]^2\\ \\
&&{\mathcal W}^{h}_{66}(\O_2)= \{ \bv_s: \bv_s|_{\O_{2,j}}\in [P_{1,1}(\O_{2,j})]^2, \,\, 
\bv_s\cdot\bnu = 0 \,\,\text{on}\,\, \G^B_2 \}\cap [C^0(\overline \O)_2]^2
\end{eqnarray*}
be the FE spaces to approximate the solid displacement, and let 
\begin{eqnarray*}
&&{\mathcal V}^h(\O)= \{\bv_f\in H(\opdiv;\O): 
\bv_f|_{\O_{j}}\in P_{1,0}(\O_j)\times P_{0,1}(\O_j),
\, \bv_f\cdot\bnu = 0 \,\,\text{on}\,\, \,\G \}\\
&&{\mathcal V}^h(\O_2)= \{\bv_f\in H(\opdiv;\O_2): 
\bv_f|_{\O_{2,j}}\in P_{1,0}(\O_{2,j})\times P_{0,1}(\O_{2,j}),
\, \bv_f\cdot\bnu = 0 \,\,\,\text{on}\,\, \G \}
\end{eqnarray*}
be the space to approximate the fluid displacement vector. 
Here $P_{s,t}$ denotes the polynomials of degree not greater than $s$ in $x_1$ 
and not greater than $t$ in $x_3$.

Then, for $(I,J) = (1, 1), (3, 3), (1, 3), (5, 5)$  let
\begin{eqnarray*}
&&{\mathcal Z}^h_{IJ}(\O) =   {\mathcal W}^h_{IJ}(\O) 
\times {\mathcal V}^h(\O).
\end{eqnarray*}
Also, let 
\begin{eqnarray*}
&&{\mathcal Z}^h_{66}(\O_2) =   {\mathcal W}^h_{IJ}(\O_2)\times {\mathcal V}^h(\O_2). 
\end{eqnarray*}

Next, for $(I,J) = (1,1), (3,3), (1,3), (5,5)$  let   
\begin{eqnarray}
\Pi^h_{IJ}: [H^{3/2}(\O)]^2
\rightarrow {\mathcal W}^{h}_{IJ}(\O)\nonumber
\end{eqnarray}
be the   interpolant operators associated with the spaces 
${\mathcal W}^{h}_{IJ}$. 
More specifically,  the degrees of freedom associated with  
$\Pi^h_{IJ} \bv$ 
are the vertexes  of the rectangles $\O_j$  and if $b$ is a common node of the  adjacent rectangles  $\O_j$ and $\O_k$ 
then  $\left(\Pi^h_{IJ} \bvarphi\right)_j(b) = \left(\Pi^h_{IJ} \bvarphi\right)_k(b) $, where  $\left(\Pi^h_{IJ} \bvarphi\right)_j$ denotes the restriction of the interpolant $\Pi^h_{IJ}\bvarphi$ of 
$\bvarphi$ to $\O_j$.

Also, let 
\begin{eqnarray}
Q^h: H_0^1(\opdiv;\O)\rightarrow {\mathcal V}^h(\O)\nonumber
\end{eqnarray}
be the projection defined by
\begin{eqnarray}
&&\left\langle(Q^h \bpsi - \bpsi)\cdot\bnu, 
1\right\rangle_{B}= 0,\quad B = \G_{jk}  \,  \, \text{or} \, 
\,  B = \G_j. \nonumber
\end{eqnarray}
The approximating properties of $\Pi^h_{IJ}$  and  $Q^h$  are
 \cite{rath,ciarlet80,nedelec80}
\begin{eqnarray}
&&\|\bvarphi- \Pi^h_{IJ} \bvarphi \|_0
+ h \|\bvarphi-\Pi^h_{IJ} \var\|_{1}
\le C  h^{3/2} \|\bvarphi\|_{3/2},\label{approx1a}\\ \nonumber\\
&&\|\bpsi- Q^h \bpsi \|_0 \le C  h\|\bpsi\|_1, 
\label{approx1c}\\ \nonumber\\
&&\|\nabla\cdot (\bpsi-  Q^h \bpsi ) \|_0 \le
C  h  \left(\|\bpsi\|_1 + \|\nabla\cdot\bpsi\|_1\right).\label{approx1d}
\end{eqnarray}
The projection operators 
\begin{eqnarray}
&&\Pi^h_{66}: [H^{3/2}(\O_2)]^2
\rightarrow {\mathcal W}^{h}_{66}(\O_2),\quad
Q^h: H_0^1(\opdiv;\O_2)\rightarrow {\mathcal V}^h(\O_2)\nonumber
\end{eqnarray}
are defined  similarly and satisfy the  approximating properties 
\eqref{approx1a}-\eqref{approx1d}.

Now, we formulate the FE procedures  to determine the stiffnesses 
$p_{IJ}$'s as follows:

\begin{itemize}
\item $p_{33}(\o)$: find $\bu^{(h,33)}\in {\mathcal Z}^h_{33}(\O)$
such that  
\begin{eqnarray}\label{fem33}
&&\Lambda(\bu^{(h,33)}, \bv) =  - \left\langle \Delta P, \bv \cdot\bnu\right\rangle_{\G^T}, 
\qquad  \forall \,\, \bv \in {\mathcal Z}^h_{33}(\O). 
\end{eqnarray}

\item  $p_{11}(\o)$: find $\bu^{(h,11)}\in {\mathcal Z}^h_{11}(\O)$
such that  
\begin{eqnarray}\label{fem11}
&&\Lambda(\bu^{(h,11)}, \bv) =  - \left\langle \Delta P, \bv \cdot\bnu\right\rangle_{\G^R}, 
\qquad  \forall \,\,  \bv \in {\mathcal Z}^h_{11}(\O).  
\end{eqnarray}


\item  $p_{13}(\o)$: find $\bu^{(h,13)}\in {\mathcal Z}^h_{13}(\O)$
such that  
\begin{eqnarray}\label{fem13}
&&\Lambda(\bu^{(h,13)}, \bv) =  - \left\langle \Delta P, \bv \cdot\bnu\right\rangle_{\G^R\cup\G^T}, 
\qquad  \forall \,\,  \bv \in {\mathcal Z}^h_{13}(\O).   
\end{eqnarray}

\item  $p_{55}(\o)$: find $\bu^{(h,55)}\in {\mathcal Z}^h_{55}(\O)$
such that  
\begin{eqnarray}\label{fem55}
&&\Lambda(\bu^{(h,55)}, \bv) =  - \left\langle {\mathbf g}, \bv_s \right\rangle_{\G\setminus \G^B}, 
\qquad  \forall \,\,  \bv \in {\mathcal Z}^h_{55}(\O).  
\end{eqnarray}

\item  $p_{66}(\o)$: find $\bu^{(h,66)}\in {\mathcal Z}^h_{66}(\O_2)$
such that  
\begin{eqnarray}\label{fem66}
&&\Lambda(\bu^{(h,66)}, \bv) =  - \left\langle {\mathbf g_2}, \bv_s \right
\rangle_{\G_2\setminus \G_2^B}, 
\qquad  \forall \,\,  \bv \in {\mathcal Z}^h_{66}(\O_2).  
\end{eqnarray}
\end{itemize}
Uniqueness for the FE  procedures \eqref{fem33}-\eqref{fem66}
can be shown  with the same argument used for the continuous case. Existence
follows from finite dimensionality.

\section{A priori error estimates}

In this section we derive  the error estimates associated with the FE procedures 
\eqref{fem33}-\eqref{fem66}.

\begin{thm}\label{thm;apriori}
Assume that  
$\bu_s^{(IJ)} \in [H^{3/2}(\O)]^2, \bu_f^{(IJ)} \in H^1(\opdiv;\O)$ 
 for $(I,J) = (1,1), (3,3), (1,3)$.  Also assume  that 
the matrix $\mathbf E$ in \eqref{def_E}  is positive definite. 
Then the following a priori error estimate holds:
\begin{eqnarray}\label{apriori_IJ}
&& \|\bu_s^{(IJ)}- \bu_s^{(h,IJ)}\|_1  + \|\bu_f^{(IJ)}- \bu_f^{(h,IJ)}\|_0
+ \|\nabla\cdot (\bu_f^{(IJ)}- \bu_f^{(h,IJ)})\|_0\\
&&\quad \le C(\o)\left[ h^{1/2} \|\bu_s^{(IJ)}\|_{3/2} 
+ h \left ( \|\bu_f^{(IJ)}\|_1 + \nabla \cdot \bu_f^{(IJ)}\|_1\right)\right].\nonumber 
\end{eqnarray}
\end{thm}
{\it Proof:} 
First we  analyze the error  for the procedure \eqref{fem33} associated with 
determination of $p_{33}$.

Recall Korn's second inequality \cite{nitsche81}:
\begin{eqnarray}\label{korn}
\sum_{l,m=1,3}\|\epsilon_{lm}(\bv)\|^2
  + \|\bv\|^2_0 \ge C_1 \|\bv\|^2_1,
\quad \forall \quad \bv\in [H^1(\O)]^2.
\end{eqnarray} 
Let $L_*$ denote the minimum eigenvalue of the matrix $\mathbf E$. 
Then use \eqref{korn} to see that the following Garding-type inequality holds:

\begin{eqnarray}\label{garding}
&&\text{Re}\left( \Lambda(\bv,\bv)\right) 
= \text{Re}\left( \left({\mathbf E} ~ \tilde e(\bv), 
\tilde e(\bv)\right)\right)\\ 
&& \qquad \ge  C_2 \left( \|\bv_s\|^2_1  + \|\nabla\cdot \bv_f\|^2_0\right)
-  C_3 \|\bv_s\|_0^2,\nonumber
\end{eqnarray}
where  
\[
C_2 = \text{min}\left(L_* C_1, 1 \right), \quad  
C_3 = L_*.
\]
Also, note that
\begin{eqnarray}\label{continuity}
&&|\Lambda(\bu, \bv) |\le C_4(\o)\left( \|\bu_s\|_1 \|\bv_s\|_1 + 
  \|\nabla\cdot \bu_f\|_0  \|\bv_s\|_1 \right. \\ 
&&  \left. \hskip2cm  + \|\bu_s\|_1 \|\nabla\cdot \bv_f\|_0 
+ \| \nabla\cdot \bu_f\|_0 \| \nabla\cdot \bv_f\|_0 +
 \|\bu_f\|_0 \|\bv_f\|_0\right)\nonumber\\
 && \hskip4cm \forall \bu=(\bu_s, \bu_f) \in {\mathcal Z}^h_{33}, \bv=(\bv_s, \bv_f) 
 \in {\mathcal Z}^h_{33}.\nonumber
\end{eqnarray}
Set
\begin{eqnarray}\label{def_e}
\be = \bu^{(33)} - \bu^{(h,33)} \equiv (\be_s, \be_f).
\end{eqnarray}
Subtract \eqref{fem33} from \eqref{weak_c33a} to obtain the error equation

\begin{eqnarray}\label{error_1}
\Lambda(\be, \bv) = 0,\quad \forall \,\,\bv\in {\mathcal Z}_{33}^h.
\end{eqnarray}
Set
\begin{eqnarray}\label{def_theta}
\bTheta^h(\bu^{(33)}) = \left(\Pi^h_{33} \bu^{(33)}_s, Q^h \bu^{(33)}_f\right)
\end{eqnarray}
and  take $\bv = \be +\bTheta^h(\bu^{(33)}) -\bu^{(33)}$ in \eqref{error_1} to obtain the equation
\begin{eqnarray}\label{error_2}
\Lambda(\be, \be) = \Lambda(\be, \bu^{(33)} - \bTheta^h(\bu^{(33)}).
\end{eqnarray}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%&&&&&&&&&&&&&&&&
Take imaginary part in  \eqref{error_2} and use the approximating properties 
\eqref{approx1a}-\eqref{approx1d} to  see that
\begin{eqnarray}\label{error_5a}
&&\o \left(\dfrac{\eta}{\kappa} \be_f, \be_f \right) = 
\text{Im}\left(\Lambda(\be, \bu^{(33)} - \bTheta^h(\bu^{(33)})\right)\\
&&\le | \Lambda(\be, \bu^{(33)} - \bTheta^h(\bu^{(33)})|\nonumber\\
&&\quad\le C_4(\o)\left[\|\be_s\|_1  \|\bu^{(33)}_s - \Pi^h_{33} \bu^{(33)}_s\|_1 
+  \| \nabla\cdot \be_f\|_0 \|\bu^{(33)}_s - \Pi^h_{33} \bu_s\|_1 \right .\nonumber\\
&&\quad  \left . +  \|\be_s\|_1 \| \nabla\cdot ( \bu^{(33)}_f - Q_h \bu^{(33)}_f))\|_0
+  \| \nabla\cdot \be_f\|_0 \| \nabla\cdot ( \bu^{(33)}_f - Q_h \bu^{(33)}_f)\|_0 \right .\nonumber\\
&&\quad  \left . + \|\be_f\|_0 \| \bu^{(33)}_f - Q_h \bu^{(33)}_f)\|_0  \right] 
\nonumber\\
&&\quad \le C_5(\o)\left[ h^{1/2} \|\bu^{(33)}_s\|_{3/2}\left( \|\be_s\|_1  + 
\|\nabla\cdot \be_f\|_0  
\right)\right.\nonumber\\
&&\quad \left . + h \|\nabla\cdot \bu^{(33)}_f\|_1 \left (\|\be_s\|_1 
+  \|\nabla\cdot \be_f\|_0 \right) 
+h \|\be_f\|_0 \|\bu^{(33)}_f\|_1  \right] \nonumber\\
&&\qquad\le \delta_1\left(\|\be_s\|_1^2 +  \|\nabla\cdot\be_f\|_0^2 \right)  
+ \delta_2 \|\be_f\|_0^2\nonumber\\ 
&&\qquad + C_{6}(\o)
\left(
h \|\bu^{(33)}_s\|^2_{3/2}  
+ h^2
\left(
\|\bu^{(33)}_f\|^2_1 + \|\nabla\cdot \bu^{(33)}_f\|^2_1
\right) 
\right) . \nonumber
\end{eqnarray}
Then, choose $\delta_2$ small enough in \eqref{error_5a} to obtain the estimate
\begin{eqnarray}\label{error_5b}
&&\|\be_f\|_0^2 \le 
 \delta_3\left(\|\be_s\|_1^2 +  \|\nabla\cdot\be_f\|_0^2 \right)  \\ 
&&\qquad + C_{7}(\o)
\left(
h \|\bu^{(33)}_s\|^2_{3/2}  
+ h^2
\left(
\|\bu^{(33)}_f\|^2_1 + \|\nabla\cdot \bu^{(33)}_f\|^2_1
\right) 
\right).\nonumber
\end{eqnarray}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Next, take real part in \eqref{error_2} and use \eqref{garding} and the argument given above 
to derive the inequality \eqref{error_5a} to obtain
\begin{eqnarray}\label{error_3}
&&C_2 \left( \|\be_s\|_1  + \|\nabla\cdot \be_f\|_0^2\right)\\ 
&&\quad \le  \text{Re} \left(\Lambda((\be, \bu - \bTheta^h)\right) 
+ C_3  \|\be_s\|^2_0 \nonumber\\
&&\quad \le  |\Lambda((\be, \bu - \bTheta^h)| 
+ C_3  \|\be_s\|^2_0 \nonumber\\
&&\quad \le C_8(\o)\left[ h^{1/2} \|\bu^{(33)}_s\|_{3/2}\left( \|\be_s\|_1 + \|\nabla\cdot \be_f\|_0   \right)\right.\nonumber\\
&&\quad \left . + h \|\nabla\cdot \bu^{(33)}_f\|_1 \left (\|\be_s\|_1 
+  \|\nabla\cdot \be_f\|_0 \right) +h \|\be_f\|_0 \|\bu^{(33)}_f\|_1  \right] 
+ C_3 \|\be_s\|^2_0\nonumber\\
&&\quad \le \epsilon \left( \|\be_s\|^2_1  +
\|\nabla\cdot \be_f\|^2_0\right) \nonumber\\
&&\quad + C_9(\epsilon,\o)\left( h \|\bu^{(33)}_s\|_{3/2} 
+ h^2\left( \|\bu^{(33)}_f\|_1^2 
+ \|\nabla\cdot \bu^{(33)}_f\|_1^2\right)\right) + 
 C_6(\epsilon,\o) \left(\|\be_s\|^2_0 + \|\be_f\|^2_0 \right)\nonumber.
\end{eqnarray}
Then, choose $\epsilon$ small enough in \eqref{error_3} to conclude that
\begin{eqnarray}\label{error_4}
&&\|\be_s\|^2_1 +  \|\nabla\cdot \be_f\|_0^2 
\le C_{10}(\epsilon,\o)
\left(
 h \|\bu^{(33)}_s\|_{3/2}  + h^2
\left( \|\bu^{(33)}_f\|_1^2 + \|\nabla\cdot \bu^{(33)}_f\|_1^2
\right)
\right)\nonumber\\
&&\qquad \qquad + C_{11}(\epsilon,\o) 
\left(\|\be_s\|^2_0 + \|\be_f\|^2_0 \right).
\end{eqnarray}
Next, we will employ a duality argument  to estimate the $\|\be\|_0$-term in the right hand side 
of \eqref{error_4}. 

Let us solve the adjoint  problem
\begin{eqnarray}\label{adjoint}
&&- i \o{\mathcal B}\bpsi - {\mathcal L^{*}}(\bpsi)= \be,\\
&&\bsig(\bpsi) \bnu \cdot \bnu = 0, \quad (x_1, x_3)\in \G^T,\label{a_bc1b}\nonumber\\
&&\bsig(\bpsi) \bnu \cdot \bchi = 0, \quad (x_1, x_3)\in \G,\label{a_bc1c}\nonumber\\
&&\bu_s\cdot\bpsi = 0, \quad (x_1, x_3)\in \G^L\cup\G^R\cup\G^B,\label{a_bc1e}\nonumber\\
&&\bu_f\cdot\bpsi=0, \quad (x_1, x_3)\in \G\label{a_bc1g}.\nonumber
\end{eqnarray}
The following regularity  will be assumed for the solution of \eqref{adjoint} \cite{kenig_88}:
\begin{eqnarray}\label{regularity}
\|\bpsi_s\|_{3/2} + \|\bpsi_f\|_1 + \|\nabla \cdot \bpsi\|_1 \le C(\o) \|\be\|_0.  
\end{eqnarray}
Thus, using integration by parts, for  any $\bv \in {\mathcal Z}^h_{33}$
\begin{eqnarray}\label{error_6}
\|\be\|_0^2 = \left(\be,- i \o B \bpsi - {\mathcal L^{*}}(\bpsi)\right)
= \Lambda(\be, \bpsi) = \Lambda(\be, \bpsi- \bv).
\end{eqnarray}
Choose $\bv = \left( \Pi^h_{33} \bpsi_s, Q^h \bpsi_f)\right)$ in \eqref{error_6},  
  use the approximating properties 
\eqref{approx1a}-\eqref{approx1d} and apply \eqref{regularity} to see that

\begin{eqnarray}\label{error_8}
&&\|\be\|_0^2 \le  C_{12}(\o) 
\left[ h^{1/2} \|\bpsi_s\|_{3/2} \left(\|\be_s\|_1 + \|\nabla\cdot\be_f\|_0\right) \right .\\ 
&&\qquad\left . + h \|\nabla\cdot \bpsi_f\|_1
\left(
\|\be_s\|_1 
+ \|\nabla \cdot \be_f\|_0  
\right) + h \|\be_f\|_0 \|\bpsi_f\|_1 \right]\nonumber\\
&&\qquad \le C_{12}(\o) \left( h^{1/2}  \left( \|\be_s\|_1 
 +  \|\nabla\cdot \be_f\|_0 \right)+ h \|\be_f\|_0\|\be\|_0\right).
\end{eqnarray}
Hence,
\begin{eqnarray}\label{error_9}
&&\|\be\|_0\le  C_{12}(\o) \left ( h^{1/2}  \left(\|\be_s\|_1 
 +  \|\nabla\cdot\be_f\|_0\right)  + h \| \be_f\|_0\right).
\end{eqnarray}
Next, use \eqref{error_5b} in \eqref{error_9} to obtain
\begin{eqnarray}\label{error_9a}
&&\|\be\|^2_0\le  C_{13}(\o) \left ( h \left(\|\be_s\|_1^2 + \|\nabla\cdot \be_f\|^2_0\right) 
 +  h^2 \delta_3 \left(\|\be_s\|^2_1 + \|\nabla \cdot \be_f\|_0^2\right)\right)\\ 
&&\qquad  + C_{14}(\epsilon,\o)
\left(
 h \|\bu^{(33)}_s\|^2_{3/2}  + h^2
\left( \|\bu^{(33)}_f\|_1^2 + \|\nabla\cdot \bu^{(33)}_f\|_1^2
\right)
\right).\nonumber
\end{eqnarray}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Employing the estimate  \eqref{error_9a} in \eqref{error_4} we conclude that  
\begin{eqnarray}\label{error_10}
&&\|\be_s\|^2_1 + \|\nabla\cdot \be_f\|_0^2 \\
&&\le C_{15}(\o)
\left(
 h \|\bu^{(33)}_s\|^2_{3/2}  + h^2
\left( \|\bu^{(33)}_f\|_1^2 + \|\nabla\cdot \bu^{(33)}_f\|_1^2
\right)
\right)\nonumber\\
&&\qquad + C_{11}(\o) \left( h \left(\|\be_s\|^2_1 
 +  \|\nabla\cdot \be_f\|^2_0 \right)+ h^2 \delta_3 \left(\|\be_s\|_1^2 
 + \|\nabla\cdot \be_f\|_0^2\right)\right).\nonumber
\end{eqnarray}
Then, choose $\delta_3$ sufficiently small in \eqref{error_10} to get
\begin{eqnarray}\label{error_11}
&&\|\be_s\|^2_1 + \|\nabla\cdot \be_f\|_0^2 \\
&&\le C_{16}(\o)
\left(
 h \|\bu^{(33)}_s\|^2_{3/2}  + h^2
\left( \|\bu^{(33)}_f\|_1^2 + \|\nabla\cdot \bu^{(33)}_f\|_1^2
\right)
\right)\nonumber\\
&&\qquad + C_{11}(\o) 
\left( h 
\left(\|\be_s\|^2_1 
 +  \|\nabla\cdot \be_f\|^2_0 
 \right)
 \right).\nonumber
\end{eqnarray}
Thus, for $h$ sufficiently small, from \eqref{error_11}  we get the inequality 
\begin{eqnarray}\label{error_12}
&&\|\be_s\|_1 +  \|\nabla\cdot \be_f\|_0 \\
&&\le C_{17}(\o)
\left(
h^{1/2} \|\bu^{(33)}_s\|_{3/2}  + h 
\left(
\|\bu^{(33)}_f\|_1 + \|\nabla\cdot \bu^{(33)}_f\|_1
\right)
\right).\nonumber
\end{eqnarray}
Finally, using \eqref{error_12} in \eqref{error_5b} we conclude 
the validity of the error estimate \eqref {apriori_IJ} 
for the solution  of \eqref{fem33}.
The error analysis for the solution of \eqref{fem11} and \eqref{fem13} 
follows with the same argument.  This completes the proof.

The error analysis for the FE procedure \eqref{fem55} and \eqref{fem66}
 requires another argument  because the solution vanishes on a set of positive measure of $\G$.
The following theorem states the estimates.

\begin{thm}\label{thm;apriori55}
Assume that  
$\bu_s^{(55)} \in [H^{3/2}(\O)]^2, \bu_f^{(55)} \in H^1(\opdiv;\O)$  
and that $\bu_s^{(66)} \in [H^{3/2}(\O_2)]^2, \bu_f^{(66)} \in H^1(\opdiv;\O_2)$.
Also, assume  that 
the matrix $\mathbf E$ in \eqref{def_E} is positive definite. 
Then for $(I,J)= (5,5), (6,6)$ the following a priori error estimate holds.
\begin{eqnarray}\label{apriori_55_66}
&& \|u_s^{(IJ)}- u_s^{(h,IJ)}\|_1  + \|u_f^{(IJ)}- u_f^{(h,IJ)}\|_0
+ \|\nabla\cdot (u_f^{(IJ)}- u_f^{(h,IJ)})\|_0\\
&&\quad \le C(\o)\left[ h^{1/2} \|u_s^{(IJ)}\|_{3/2} 
+ h \left ( \|u_f^{(IJ)}\|_1 + \nabla \cdot u_f^{(IJ)}\|_1\right)\right].\nonumber 
\end{eqnarray}
\end{thm}
{\it Proof:}
Set
\begin{eqnarray}\label{def_e55}
\be = \bu^{(55)} - \bu^{(h,55)} = (\be_s, \be_f).
\end{eqnarray}
Subtract \eqref{fem55} from \eqref{weak_c55} to obtain the error equation

\begin{eqnarray}\label{error_155}
\Lambda(\be, \bv) = 0,\quad \forall \bv\in {\mathcal Z}_{55}^h .
\end{eqnarray}
Set
\begin{eqnarray}\label{def_theta55}
\bTheta^h(\bu^{(55)}) = \left(\Pi^h_{55} \bu^{(55)}_s, Q^h \bu^{(55)}_f\right)
\end{eqnarray}
and  take $\bv = \be +\bTheta^h(\bu^{(55)}) -\bu^{(55)}$ in \eqref{error_155}
to obtain 
\begin{eqnarray}\label{error_255}
\Lambda(\be, \be) = \Lambda(\be, \bu^{(55)} - \bTheta^h(\bu^{(55)}).
\end{eqnarray}
Next, recall that 
\begin{eqnarray}\label{equiv_norm1}
\||\bv\|| = \left( \sum_{l,m=1,3}\|\epsilon_{lm}(\bv)\|^2 \right)^{1/2} 
\end{eqnarray}
defines a norm on  $\{\mathcal W^h_{55}$  equivalent to the $H^1$-norm \cite{ciarlet80},  i.e., 
there exist  positive  constants $C_{11}, C_{12}$ such that
\begin{eqnarray}\label{equiv_norm2}
C_{18} \|\bv\|_1\le \||\bv \|| \le C_{19} \|\bv\|_1, \quad \forall \,\, \bv\in {\mathcal W}^h_{55} .
\end{eqnarray}
Hence, 
\begin{eqnarray}\label{error_355}
&&\text{Re}\left( \Lambda(\be,\be)\right) 
= \text{Re}\left( \left({\mathbf E} ~ \tilde e(\be), 
\tilde e(\be)\right)\right)\\ 
&& \qquad \ge L_* \left( \||\be_s\||^2 + \|\nabla\cdot \be_f\|^2_0 \right)
\ge \nonumber L_*\left(C_{18} \|\be_s\|_1^2  + \|\nabla\cdot \be_f\|^2_0 \right) .
\end{eqnarray}
Then, take real part in \eqref{error_255} and use the approximating properties 
\eqref{approx1a}-\eqref{approx1d} to obtain
\begin{eqnarray}\label{error_455}
&&L_*\left(C_{18} \|\be_s\|_1^2  + \|\nabla\cdot \be_f\|^2_0 \right)\\ 
&&\qquad \le \text{Re}\left( \Lambda(\be, \bu^{(55)} - \bTheta^h(\bu^{(55)}) \right)
\le | \Lambda(\be, \bu^{(55)} - \bTheta^h(\bu^{(55)})|\nonumber\\
&&\qquad\le C_{20}(\o)\left[ h^{1/2} \|\bu^{(55)}_s\|_{3/2}
\left( \|\be_s\|_1  + \|\nabla\cdot \be_f\|_0   \right)\right.\nonumber\\
&&\quad \left . + h \|\nabla\cdot \bu^{(55)}_f\|_1 \left (\|\be_s\|_1 +  
\|\nabla\cdot \be_f\|_0 \right) +h \|\be_f\|_0 \|\bu^{(55)}_f\|_1  \right].\nonumber\\
&&\qquad \le \epsilon\left(\|\be_s\|_1^2 + \|\be_f\|^2_0 
+ \|\nabla\cdot \be_f\|^2_0\right)\nonumber\\ 
&&\qquad + C_{21}(\o)
\left(
h \|\bu^{(55)}_s\|^2_{3/2}  
+ h^2
\left(
\|\bu^{(55)}_f\|^2_1 + \|\nabla\cdot \bu^{(55)}_f\|^2_1
\right) 
\right) . \nonumber
\end{eqnarray}
Next, taking  imaginary part in  \eqref{error_255} and  repeating the argument leading to \eqref{error_5b}
we see that the $\|\be_f\|_0$-term in the right-hand side of \eqref{error_455} satisfies the estimate
\begin{eqnarray}\label{error_655}
&&\|\be_f\|^2_0 \le  \delta_3\left(\|\be_s\|_1^2 +   \|\nabla\cdot \be_f\|^2_0\right)\\ 
&&\qquad + C_{7}(\o)\left(
h \|\bu^{(55)}_s\|^2_{3/2}  
+ h^2
\left(
\|\bu^{(55)}_f\|^2_1 + \|\nabla\cdot \bu^{(55)}_f\|^2_1
\right) 
\right) . \nonumber
\end{eqnarray}
Next  using \eqref{error_655} in \eqref{error_455}, we obtain
\begin{eqnarray}\label{error_755}
&&L_*\left(C_{18} \|\be_s\|_1^2  + \|\nabla\cdot \be_f\|^2_0 \right)\\
&&\qquad \le \epsilon\left(\|\be_s\|_1^2  
+ \|\nabla\cdot \be_f\|^2_0\right)\nonumber\\ 
&&\qquad + C_{21}(\o)\left( h \|\bu^{(55)}_s\|^2_{3/2}  
+ h^2\left(\|\bu^{(55)}_f\|^2_1 + \|\nabla\cdot \bu^{(55)}_f\|^2_1\right) 
\right)\nonumber\\
&&\qquad + \epsilon \left(\delta_3\left(\|\be_s\|_1^2 +  
\|\nabla\cdot \be_f\|^2_0\right) \right.\nonumber\\ 
&&\left .\qquad + C_{7}(\o)\left(
h \|\bu^{(55)}_s\|^2_{3/2}  
+ h^2
\left(
\|\bu^{(55)}_f\|^2_1 + \|\nabla\cdot \bu^{(55)}_f\|^2_1
\right) 
\right)\nonumber    \right) 
\end{eqnarray}
Thus,  take $\epsilon$ and $\delta_3$ small enough in \eqref{error_755} to derive the inequality
\begin{eqnarray}\label{error_855}
&&\|\be_s\|_1^2  + \|\nabla\cdot \be_f\|^2_0 \le C_{22}(\o)
\left [ h\|\bu_s^{(55)}\|^2_{3/2} \right .\\
&&\left .\hskip5cm + h^2\left(\|\bu^{(55)}_f\|^2_1 
+  \|\nabla\cdot \bu^{(55)}_f\|^2_1 \right)\right].\nonumber
\end{eqnarray}
Finally using \eqref{error_855} in \eqref{error_655} we conclude the validity of the error
estimate \eqref{apriori_55_66} for $(I,J)= (5,5)$. The proof for $(I,J)= (6,6)$ is identical.
This completes the proof.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Numerical experiments}

The FE procedures described above are implemented  to determine 
the five complex stiffnesses 
$p_{IJ}(\o)$ as a function  of frequency  
and the corresponding  phase velocities and  dissipation coefficients 
as indicated in Appendix B. 
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 moduli of the samples are determined by 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}
Porosity is $\phi = 0.25$ in the background and $\phi = 0.5$ in the fractures. Using 
\eqref{krief} we obtain   $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, giving 
$\kappa = 0.247$ D  in the background and $\kappa = 4.44$ D in the fractures.

The first numerical experiments show the validation of the FE 
 procedure by comparison  with the  analytical solution given in \cite{KM_11}
 which is included in  Appendix A. 
 In this example, we consider a brine saturated sample, with  brine having a density of 1040 kg/m$^3$, a viscosity of
 0.0018 Pa s and a bulk modulus of 2.25 GPa.
 
 Figures \ref{fig1_vp11_vp33_freq} and \ref{fig2_qp11_qp33_freq} show the compressional wave phase velocity 
 and dissipation factor, respectively,  as a function of frequency,  
 in the direction  parallel (squares and solid lines) and normal (circles 
 and solid lines) to the fractures. 
The solid lines indicate the theoretical values, while symbols 
indicate the FE solution. 
It can be observed a perfect fit  of the FE solution to the  analytical values in the whole frequency range displayed.
\vskip1cm
$$
< < \text{Figure 1  }> >
$$
\vskip1cm 


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

The phase  velocities of the qP, qSV and SH waves and dissipation factors of the qP and qSV 
waves as a function of the propagation angle are represented
in Figures \ref{fig3_qP_qSV_SH_veloc_angle} and \ref{fig4_qP_qSV_SH_aten_angle}, respectively, where the frequency is 300 Hz. 
Symbols and solid lines indicate 
the FE and theoretical values, respectively. Again an excellent 
match between the numerical and  analytical curves is obtained for all angles.
Here and in the following figures, the propagation angle is 
understood to be given 
with respect to the symmetry axis, so that 0 degrees and 90  degrees correspond 
to waves arriving normal and parallel to the fracture layering, respectively. 
Figures \ref{fig3_qP_qSV_SH_veloc_angle} and \ref{fig4_qP_qSV_SH_aten_angle} show that anisotropy induced  by fractures is noticeable for the phase velocities of all waves. 
On the other hand, the qP curves show strong attenuation 
for waves arriving normal to the 
fracture layering. The qSV  wave has no loss along the
directions parallel and perpendicular to the layering plane, showing maximum 
attenuation at about 45 degrees.

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


$$
< < \text{Figure 4  }>>
$$
\vskip1cm
  
Next, we  present two experiments for which no analytical solutions are available.   
 The first experiment considers   the case  of a patchy brine-CO$_2$ saturated fractured sample, 
 with  CO$_2$ having density  500 kg/m$^3$, viscosity 
 2.10$^{-5}$ Pa s and bulk modulus  22.5 MPa.
To generate a patchy CO$_2$-brine saturation we use the von Karman self-similar 
correlation function  for which the spectral density is given
by \cite{Frankel}
\begin{equation}\label{vonKarman}
S_d(r_x,r_z)= N_0 (1 + R^2 a ^2)^{-(H + E/2)}. 
\end{equation}
Here, $R=\sqrt{r_x^2+r_z^2}$ is the radial wavenumber, $a$ the correlation length, 
$H$ is a self-similarity coefficient ($0 < H < 1$), $N_0$ is a normalization constant 
and $E$ is the euclidean dimension.
The von Karman correlation  \eqref{vonKarman}
 describes a self-affine, fractal processes of fractal dimension
$D=E+1-H$ at a scale smaller than $a$. We choose $E$ = 2,  $D$ = 2.2 
and $a$ to be  2.5 \% of the domain size.  Once a  continuous fractal 
distribution of brine is obtained over  the 160 $\times$ 160 mesh,   
by properly chosen threshold values $S_b^*$, for each 
cell with brine saturation below  and above $S_b^*$ 
we assign to that cell  either full   CO$_2$ or full  brine saturation, respectively.
In this way 
we generated two different 
patchy CO$_2$-brine distributions of overall  10 $\%$ and 30 $\%$  CO$_2$ saturation. 
Figure \ref{fig5_patchy_10_CO2_brine} displays the obtained  patchy CO$_2$-brine saturation at 
10  \%  CO$_2$ saturation. 

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


Figures \ref{fig6_qP_veloc_brine_10_30_CO2_angle} and  \ref{fig7_qsV_veloc_brine_10_30_CO2_angle}   show the   phase velocities of the qP and qSV waves at 300 Hz for full brine saturation, 
10 $\%$ and 30 $\%$ patchy  CO$_2$-brine  saturation  as function of the 
propagation angle, while Figures \ref{fig8_qP_aten_brine_10_30_CO2_angle} and  \ref{fig9_qSV_aten_brine_10_30_CO2_angle} 
show the corresponding dissipation factors.  

It can be seen that for qP waves, patchy saturation enhances  the 
strong velocity anisotropy  caused by fractures observed in the brine saturated case, with values 
 depending on the overall CO$_2$ 
saturation.   Maximum differences  with respect to the full brine saturation case are
 at 10 $\%$ patchy saturation 
and angles normal to the fracture layering. 

On the other hand,  the qSV phase velocity anisotropy is less affected by patchy saturation. 
For both the qP and qSV waves 
we observe lower velocities when patches are present, with decreasing values for  
increasing CO$_2$ overall saturation. SH waves phase velocities are not affected by the 
presence of  CO$_2$ patches, so the corresponding curve is identical to that  
for the full brine saturation in Figure \ref{fig3_qP_qSV_SH_veloc_angle} and it is omitted.

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

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

Concerning the dissipation factors, the patchy saturation enhances 
attenuation anisotropy of the qP waves for all the angles and it is strong at 10 \% patchy saturation 
for waves arriving normal to the fracture layering plane and up to 30 degrees. Above 60 degrees, for the two cases of patchy saturation qP waves show an almost  constant dissipation factor independently 
 of the value of CO$_2$ saturation. This is opposite to the case of full brine saturation, for which dissipation vanishes for angles above 60 degrees.

For qSV waves, attenuation  anisotropy is null for waves arriving normal or parallel to the fracture layering plane and 
strong for angles between 30 and  
60 degrees. Also, when compared with full brine saturation,  dissipation factors are almost unaffected for 10 \% patchy saturation but noticeably reduced for the  30 \%  CO$_2$ case.

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

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

Figures \ref{fig10_pressure_33_10_CO2} and \ref{fig11_pressure_33_30_CO2} show the fluid pressure distribution for compressions normal to the 
fracture layering ($p_{33}$ experiment) 
at frequency 300 Hz and  for  10 \% and 30 \%  patchy CO$_2$ saturation, respectively.
 The figure for 10 \%  CO$_2$ saturation  has larger regions of high pressure gradients  
than  the 30 \%  one , which explains the lower values of the dissipation factors for qP and qSV waves  for 30 \%  CO$_2$ saturation (squares) as compared with the corresponding 
curves for 10 \%  CO$_2$ saturation (circles) in Figures \ref{fig8_qP_aten_brine_10_30_CO2_angle} and  \ref{fig9_qSV_aten_brine_10_30_CO2_angle}. 
\vskip1cm
$$
< < \text{Figure 10  }>>
$$
\vskip1cm

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

The second experiment considers the case of the  brine saturated sample with uniform 
background as in the 
first set of experiments  but with 
 fractures having  binary fractal variations in their petrophysical 
properties, which is another case not described by the theory in \cite{KM_11}.
To generate such sample, first we generated a continuous fractal distribution of porosity 
over the mesh by using the spectral density in \eqref{vonKarman}  with average 
$\overline \phi$ equal to 
 0.5,  fractal dimension $D$ = 2.3 and correlation length $a =$  0.3  in a scale of 10. 
In this way, we obtain minimum and maximum values of porosity of   0.32  and 0.646, respectively.
Permeability is then computed using \eqref{perm}, giving values between 0.7  Darcy 
and 19 Darcy, with an average of 4.9 Darcy.
 The binary fractals associated with porosity and permeability 
are then computed as follows.  For each cell in the computational mesh, 
if the value of porosity is smaller (respectively, bigger) than the average $\overline \phi$, 
we assign to that cell the minimum (respectively, maximum) porosity value. 
The procedure is repeated to obtain a binary fractal permeability distribution. 
Binary fractal distributions of bulk modulus $K_m$ and shear modulus 
$\mu$ are determined by using \eqref{krief} and the computed binary fractal porosity field. 

Taking 1D restrictions of these four binary fractal distributions for porosity $\phi$, permeability $\kappa$, and bulk and shear moduli $K_m$ and $\mu$, we obtain the binary fractal petrophysical properties defining the  brine saturated fractures.

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

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

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

Figures \ref{fig12_veloc_qP_uniform_fractal_fractures},  
\ref{fig13_veloc_qP_uniform_fractal_fractures} and 
\ref{fig14_veloc_SH_uniform_fractal_fractures} display the phase velocities of the qP qSV and SH waves at 300 Hz as a function of the propagation angle for a brine saturated sample with uniform and binary fractal variations in the petrophysical properties of the fractures. Figures \ref{fig15_aten_qP_uniform_fractal_fractures}, 
\ref{fig16_aten_qSV_uniform_fractal_fractures} and 
\ref{fig17_aten_SH_uniform_fractal_fractures} show the corresponding dissipation factors.
It can be noticed a strong increase in  velocity anisotropy for the three waves with respect to the case of uniform fractures. 
In particular, the qSV wave 
suffers a velocity reduction for almost all angles except near 40 degrees and the phase velocity
of the SH wave decreases by about  75 \%. 
Also, qP waves have a velocity increase for waves arriving at angles close to the normal to the fracture layering, and  velocity decreases for all other angles.
Regarding the dissipation factors, the qP waves show an increase in attenuation for angles up to 50 degrees, while the qSV waves also show increasing attenuation for angles above 45 degrees. 
The most significant changes are observed in the SH dissipation factors, which change  from vanishing  dissipation  to a continuous increase in dissipation for increasing angles.

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

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

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

\section{Conclusions}

We have presented a set of finite element harmonic experiments to 
determine the five complex and frequency-dependent 
stiffnesses of the TIV medium equivalent to 
a fractured fluid-saturated porous material, where the fractures 
are modeled as very thin, highly permeable and compliant porous layers.
The procedure allow us to compute the wave
velocities and quality factors at the macroscale 
as a function of frequency and propagation angle. 
The proposed  methodology  is based on the 
solution of the diffusive Biot equations in the space-frequency domain
to simulate harmonic compressibility and shear tests.
The methodology is validated  against a  theory valid 
at long wavelengths  for homogeneous layers and 
fluid flow normal to the fracture layering. 

Then, the experiments are applied for the cases of patchy brine-CO$_2$ saturation and 
a brine saturated sample of uniform background and fractures with fractal variations in their  
petrophysical properties. 
For the case of a patchy brine-CO$_2$ saturated sample, the  experiments show that 
both for qP and qSV waves 
the presence of patches  
enhances the strong seismic velocity and attenuation 
anisotropy induced by the fractures. 
Also,  SH waves show moderate velocity anisotropy, whose values are not 
affected by the presence of the patches of CO$_2$. 
In the last experiment, considering a brine saturated sample of uniform background
and fractures of having fractal variations in petrophysical properties, it is observed a noticeable 
increase in velocity and Q anisotropy for all the waves, with higher dissipation factors 
as compared with  the case of  fractures of uniform properties. In particular, the SH dissipation factor 
shows positive and increasing values with increasing angles, as opposite to the uniform fractures case, which exhibits 
vanishing attenuation at all angles.
   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Aknowledgements}
 The work of J. E. 
Santos was partially funded by the GRANT PIP 112-20080100952 from CONICET, 
Argentina. J. M. Carcione has been partially funded by the CO2Monitor project. 
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\newpage
\begin{thebibliography}{999}
%%%%%%%%%%%%  1 %%%%%%%%%%%%%%
\bibitem{GU_09}
B. Gurevich, M.  Brajanovski, R. J. Galvin, T. M. M\"uller and
J. Toms-Stewart,  P-wave dispersion and attenuation
in fractured and porous reservoirs--poroelasticity approach,
Geophys. Prospect.,  57 (2009) 225-237.

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

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

%%%%%%%%%%%%%%%%%%%   4 

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


%%%%%%%%%%%%%%  5
\bibitem{CP_06}
J. M. Carcione  and S. Picotti, P-wave seismic attenuation by slow-wave diffusion,
Effects of inhomogeneous rock properties, Geophysics 71 (2006)  O1-O8. 

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

%%%%%%%%%%%%%%%   7

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

%%%%%%%%%%%%%%%   8
\bibitem{CGSP_13}
J. M. Carcione, B. Gurevich, J. E. Santos   and S. Picotti, 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{SRR_09}
J. E. Santos, J. G.  Rubino and C. L.  Ravazzoli, 
A numerical upscaling procedure to estimate  effective bulk and shear
moduli in heterogeneous fluid-saturated porous media,
Comput. Methods Appl. Mech. Engrg. 198 (2009) 2067-2077.

%%%%%%%%%%%%%%%%%%   10

\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 (2010) L06302.

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

\bibitem{SPC_12} 
J. E. Santos, S. Picotti and Carcione, 
Evaluation of the stiffness tensor of  a fractured medium with harmonic experiments, 
Comput. Methods Appl. Mech. Engrg. 247-248 (2012) 130-145. 


%%%%%%%%%%%%%   12
\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 (2007) SM293-SM300.

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

%%%%%%%%%%%%%%%%%%   14
 

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


%%%%%%%%%%%   15 
\bibitem{carcione_07}
J. M. Carcione, 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, 2007.


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

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

\bibitem{kolsky63}
Kolsky, H., 1963. {\it Stress waves in solids},
Dover publications, New York.

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

\bibitem{adams_75} 
R. A. Adams,  Sobolev Spaces, Academic Press, 1975.

%%%%%%%%%%%%%%%  19

\bibitem{kenig_88}
B. E. J. Dahlberg, C. E. Kenig and G. C. Verchotta,
Boundary value problems for the systems of elastostatics in
Lipschitz domains,
Duke Math. J., 57 (1988) 795-818.

%%%%%%%%%%%%%%%  19

\bibitem{ciarlet80}
P. G. Ciarlet, The Finite Element Method for Elliptic Problems, 
Noth Holland, 1980.

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

\bibitem{nitsche81}
J. A. Nitsche, On Korn's second inequality,
RAIRO, 15 (1981) 237-249. 

%%%%%%%%%%%%%%%  21


\bibitem{nedelec80}
J. C. Nedelec,
Mixed finite elements in $R^3$, Numer. Math. 35 (1980) 315-341.


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

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


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

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

%%%%%%%%%%%%%%   24
 
 \bibitem{Frankel}
A.~Frankel, R.~W. Clayton, Finite difference simulation of seismic wave
  scattering: implications for the propagation of short period seismic waves in
  the crust and models of crustal heterogeneity, J. Geophys. Res. 91 (1986)
  6465 -- 6489.


\end{thebibliography}

\newpage 

\listoffigures

\newpage

\appendix


\section{Mesoscopic-flow attenuation theory for anisotropic poroelastic media}
White's mesoscopic attenuation theory of interlayer flow \cite{white_75}, \cite{CP_06}
 describes
the equivalent viscoelastic medium of a stack of 
two thin alternating porous layers of thickness $d_1$ 
and $d_2$, such that the period of the stratification is $d = d_1 + d_2$. The 
theory gives the complex and frequency dependent stiffness $p_{33}$. 
White model 
has been generalized in \cite{KM_11} by Krzikalla and M\"uller
to anisotropic media, i.e., 
they have obtained the five 
stiffnesses of the equivalent transversely isotropic medium, 
denoted by $p_{IJ}$. The stress-strain relations is 
given by equations (\ref{st_GL_1})-(\ref{st_GL_6}) and  
\begin{equation}\label{stre-stra}
p_{IJ}( \omega )  = c_{IJ} + \left( \frac{ c_{IJ} - c_{IJ}^r}{c_{33} - c_{33}^r} \right)  [p_{33} (\omega)  - c_{33}] ,
\end{equation}
where $c_{IJ}^r$ and $c_{IJ}$ are the relaxed and unrelaxed stiffnesses. 

According to Gelinsky and Shapiro \cite{GS_97} [their eq. (14)], 
the quasistatic or relaxed effective constants of a stack of poroelastic layers are
\begin{equation}\label{eff}
\begin{array}{l}
c_{66}^r = B_1^\ast = \langle \mu \rangle , \\
c_{11}^r - 2 c_{66}^r = c_{12}^r  = B_2^\ast = 2 \left \langle \myfrac{\lambda_m \mu}{E_m} \right \rangle + 
\left \langle \myfrac{\lambda_m }{E_m} \right \rangle^2  \left \langle \myfrac{1}{E_m} \right \rangle^{-1} \\
\hspace{3.6cm} +  \myfrac{{B_6^\ast}^2}{B_8^\ast}, \\
c_{13}^r  = B_3^\ast = \left \langle \myfrac{\lambda_m}{E_m} \right \rangle   
\left \langle \myfrac{1}{E_m} \right \rangle^{-1} +  \myfrac{B_6^\ast B_7^\ast}{B_8^\ast}, \\
c_{33}^r  = B_4^\ast =   
\left \langle \myfrac{1}{E_m} \right \rangle^{-1} +  \myfrac{{B_7^\ast}^2}{B_8^\ast} \\ 
%Norris
\hspace{1.3cm}
= \left[ \left \langle \myfrac{1}{E_m} \right \rangle - \left \langle \myfrac{\alpha}{E_m} \right \rangle^2  
\left \langle \myfrac{E_G}{M E_m} \right \rangle^{-1} \right]^{-1} , \\
c_{55}^r = B_5^\ast = \langle \mu^{-1} \rangle^{-1} , \\
\ \ \ \ \ B_6^\ast =  - B_8^\ast \left( 
2 \left \langle \myfrac{\alpha \mu}{E_m} \right \rangle + \left \langle \myfrac{\alpha}{E_m} \right \rangle
\left \langle \myfrac{\lambda_m}{E_m} \right \rangle \left \langle \myfrac{1}{E_m} \right \rangle^{-1} \right)  , \\
\ \ \ \ \ B_7^\ast =  - B_8^\ast \left \langle \myfrac{\alpha}{E_m} \right \rangle \left \langle \myfrac{1}{E_m} \right \rangle^{-1}   , \\
\ \ \ \ \ B_8^\ast =
\left[ \left \langle \myfrac{1}{M} \right \rangle + \left \langle \myfrac{\alpha^2}{E_m} \right \rangle - 
\left \langle \myfrac{\alpha}{E_m} \right \rangle^2 \left \langle \myfrac{1}{E_m} \right \rangle^{-1} 
\right]^{-1} ,
\end{array}
\end{equation}
where 
\begin{equation}\label{dry}
\begin{array}{l}
\lambda_m = K_m - \myfrac{2}{3} \mu  \ \ \ {\rm and} \ \ \
E_m = K_m + \myfrac{4}{3} \mu 
\end{array}
\end{equation}
and we have also reported the notation of that paper for clarity.
In the case of no interlayer flow, i.e., the unrelaxed regime, the stiffnesses are 
\begin{equation}\label{effnf}
\begin{array}{l}
c_{66} =  c_{66}^r , \\
c_{11} - 2 c_{66} = c_{12} = 2 \left \langle \myfrac{(E_G - 2 \mu)  \mu}{E_G} \right \rangle + 
\left \langle \myfrac{E_G - 2 \mu}{E_G} \right \rangle^2  \left \langle \myfrac{1}{E_G} \right \rangle^{-1}, \\
c_{13} = \left \langle \myfrac{E_G - 2 \mu}{E_G} \right \rangle   
\left \langle \myfrac{1}{E_G} \right \rangle^{-1} , \\
c_{33} = 
\left \langle \myfrac{1}{E_G} \right \rangle^{-1} , \\
c_{55} =  c_{55}^r 
\end{array}
\end{equation}
[Gelinsky and Shapiro (1997), eq. (15)], where  
\begin{equation}\label{wet}
E_G = E_m + \alpha^2 M, 
\end{equation} 
and $M$ is given in (\ref{mod.3}).

Finally, the P-wave modulus $p_{33}$ is \cite{white_75}, also see 
 in \cite{carcione_07}
\begin{equation}\label{p33}
p_{33} = \left[ \frac{1}{c_{33}} 
+ \frac{2 (r_2 - r_1)^2} {\ii \omega (d_1 + d_2) (I_1 + I_2)} \right]^{-1},
\end{equation}
where
\begin{equation}\label{rn}
r = \frac{\alpha M}{E_G}
\end{equation} 
and
\begin{equation}\label{In}
I =  \frac{\eta}{\kappa a} {\rm coth} 
\left( \frac{a d}{2} \right) , \ \ \ a = \sqrt{\frac{ \ii \omega \eta E_G}{\kappa M E_m}} , 
\end{equation}
for each single layer.

The main assumption in \cite{KM_11} is  that  
 the fluid-flow direction is perpendicular to the fracture 
 layering 
and that the relaxation behavior is described by a single
relaxation function or stiffness, i.e., $p_{33} (\omega)$. 
Thus the theory is valid for 
plane layers and  can not be used when 
 2D or 3D heterogeneities are present.


\section{Wave velocities and quality factors}

We consider homogeneous viscoelastic waves \cite{carcione_07}.
The complex velocities are the key quantity to obtain the
wave velocities and quality factor of the equivalent anisotropic medium.
They are given by 
\begin{equation}\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} , \\
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{array} 
\end{equation}
where $\bar \rho$ is the average 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 symmetry axis,
and the three velocities correspond to the qP, qS and SH waves, respectively.
The phase velocity is given by
\begin{equation} \label{vph}
v_p = \left[{\rm Re} \left( \frac { 1 }{ v } \right) \right]^{-1},
\end{equation}
while the  quality factor is given by
\begin{equation} \label{qfacvis1}
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 of the qP and qSV waves is given by
%\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}
%[eq. 6.158] in \cite{carcione_07}, where
%\begin{equation} \label{40an}
%\tan{\psi} = \frac {{\rm Re} (\beta^\ast X + \xi^\ast W)}
%{{\rm Re} (\beta^\ast W + \xi^\ast Z)} ,
%\end{equation}
%defines the angle between the energy-velocity vector and the $z$-axis,
%\begin{equation} \label{19an}
%\begin{array}{l}
%\beta = \sqrt{ A \pm B} , \\
%\xi = \pm {\rm pv} \sqrt{ A \mp B} , \\
%B = p_{11} l_1^2 - p_{33} l_3^2 + p_{55} \cos 2 \theta ,
%\end{array}
%\end{equation}
%where the upper and lower signs correspond 
%to the qP and qS waves, respectively. Moreover,
%\begin{equation} \label{23an}
%\begin{array}{l}
%W = p_{55} ( \xi l_1 + \beta l_3 ) , \\
%X = \beta p_{11} l_1 + \xi p_{13} l_3 , \\
%Z =  \beta p_{13} l_1 + \xi p_{33} l_3
%\end{array}
%\end{equation}
%where ``pv" denotes the principal value, which has
%to chosen according to established criteria.

%On the other hand, the energy velocity of the SH wave is
%\begin{equation}\label{evelSHav}
%{\bf v}_e =
%\frac{1}{\bar \rho v_p} \left( l_1 c_{66} \hat {\bf e}_1+
%l_3 c_{55} \hat {\bf e}_3 \right)
%\end{equation}
%and
%\begin{equation} \label{40anSH}
%\tan{\psi} = \left( \frac { c_{66} }{ c_{55}} \right) \tan \theta
%\end{equation}
%(Carcione 2007; eq. 1.148), since $p_{55}$ and $p_{66}$ are real quantities. 

%In general, the phase velocity is related to the energy velocity by
%\begin{equation} \label{58an}
%v_p = v_e \cos (\psi - \theta) ,
%\end{equation}
%where $v_e = |{\bf v}_e|$.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


\newpage

%\end{document}

\setcounter{figure}{0}
\renewcommand{\thefigure}{\arabic{figure}}



\begin{figure} [!h]
\begin{center}
\includegraphics[width=16cm]
{fig1_vp11_33_fe_and_analytical_brine_with_fractures.eps}

\vskip1cm

\caption{P-wave velocity as a function of frequency in the direction parallel (squares and solid lines) and normal (circles and solid lines) 
to the fracture plane. The solid lines
indicate the theoretical values.}
\label{fig1_vp11_vp33_freq} 
\end{center}
\end{figure}

\newpage

\begin{figure}[!h]
\begin{center}
\includegraphics[width=16cm]
{fig2_1000_over_q11_q33_fe_and_analytical_brine_with_fractures.eps}

\vspace{1cm}

\caption{ P-wave dissipation factors  as a function of frequency in the 
direction parallel (squares and solid lines) and normal (circles and solid lines) 
to the fracture plane. The solid lines indicate the theoretical values.}
\label{fig2_qp11_qp33_freq}
\end{center}
\end{figure}

\newpage

\begin{figure} [!h]
\begin{center}
\includegraphics[width=16cm]
{fig3_phase_velocities_qP_qSV_SH_FE_and_analytical_300Hz_versus_angle.eps}

\vspace{1cm}

\caption{Phase velocities for qP (circles), qSV (squares) and SH (triangles up) 
waves as function of the propagation angle. Frequency is 300 Hz.  
The solid lines indicate the theoretical values.}
\label{fig3_qP_qSV_SH_veloc_angle}
\end{center}
\end{figure}

\newpage

\begin{figure}[!h]
\begin{center}
\includegraphics[width=16cm]
{fig4_1000_over_qP_qSV_fe_and_analytical_300Hz_angle.eps}


\vskip1cm


\caption{Dissipation factors  for qP (circles) and  qSV (squares)  
waves as function of the propagation angle. Frequency is 300 Hz.
The solid lines indicate the theoretical values.}
\label{fig4_qP_qSV_SH_aten_angle}
\end{center}
\end{figure}




\newpage

\begin{figure}[!h]
\begin{center}
\includegraphics[width=16cm]
{fig5_matnum_patchy_gnu.dat_10_CO2_grey.eps}

\vskip3cm
\caption{Patchy CO$_2$-brine saturation. White zones correspond to full 
CO$_2$ saturation, and  black ones to  full brine saturation. 
The overall CO$_2$ saturation is 10 \%.}
\label{fig5_patchy_10_CO2_brine}
\end{center}
\end{figure}




\newpage
\begin{figure}[!h]
\begin{center}
\includegraphics[width=16cm]
{fig6_phase_velocitiy_qP_waves_brine_10_and_30_percent_patchy_CO2_saturation.eps}

\vskip3cm
\caption{qP phase velocities for brine, 10 \% and 30 \% patchy CO$_2$-brine saturation as function of the propagation angle. Frequency is 300 Hz.}
\label{fig6_qP_veloc_brine_10_30_CO2_angle}
\end{center}
\end{figure}


\newpage



\begin{figure}[!h]
\begin{center}
\includegraphics[width=16cm]
{fig7_phase_velocity_qSV_waves_brine_10_and_30_percent_patchy_saturation_300Hz_versus_angle.eps}

\vspace{1cm}

\caption{qSV  phase velocities for brine, 10 \%  and 30 \%  percent patchy  CO$_2$-brine  saturation
as function of the propagation angle. Frequency is 300 Hz.}
\label{fig7_qsV_veloc_brine_10_30_CO2_angle}
\end{center}
\end{figure}


\newpage



\begin{figure}[!h]
\begin{center}
\includegraphics[width=16cm]
{fig8_1000_over_qp_brine_10_patchy_CO2_30_patchy_CO2.eps}

\vspace{1cm}


\caption{Dissipation factors for brine, 10 \%  and 30 \%  patchy  CO$_2$-brine saturation 
as function of the propagation angle. Frequency is 300 Hz.}
\label{fig8_qP_aten_brine_10_30_CO2_angle}
\end{center}
\end{figure}

\newpage

\begin{figure}[!h]
\begin{center}
\includegraphics[width=16cm]
{fig9_1000_over_qsv_brine_10_patchy_CO2_30_patchy_CO2.eps}

\vspace{1cm}

\caption{Dissipation factors of qSV waves for brine, 10 \% and 30 \% 
 patchy  CO$_2$-brine  saturation 
as function of the propagation angle. Frequency is 300 Hz.}
\label{fig9_qSV_aten_brine_10_30_CO2_angle}
\end{center}
\end{figure}


\newpage

\begin{figure}[!h]
\begin{center}
\includegraphics[width=16cm]
{fig10_pressure_normal_to_fracture_plane_10_percent_CO2_grey.eps}

\vspace{1cm}

\caption{Fluid pressure for normal compression to the fracture plane at  10 \%  
patchy CO$_2$-brine saturation. Frequency is 300 Hz.}
\label{fig10_pressure_33_10_CO2}
\end{center}
\end{figure}

\newpage

\begin{figure}[!h]
\begin{center}
\includegraphics[width=16cm]
{fig11_pressure_normal_to_fracture_plane_30_percent_CO2_grey.eps}

\vspace{1cm}

\caption{Fluid pressure for normal compression to the fracture plane at  30 \% 
patchy  CO$_2$-brine saturation. Frequency is 300 Hz.}
\label{fig11_pressure_33_30_CO2}
\end{center}
\end{figure}


\newpage

\begin{figure}[!h]
\begin{center}
\includegraphics[width=16cm]
{fig12_phase_velocity_qP_waves_uniform_and_fractal_fractures_300Hz_angle.eps}

\vspace{1cm}

\caption{Phase velocities of qP waves for a brine saturated sample with 
uniform and binary fractal fractures.}
\label{fig12_veloc_qP_uniform_fractal_fractures}
\end{center}
\end{figure}

\newpage

\begin{figure}[!h]
\begin{center}
\includegraphics[width=16cm]
{fig13_phase_velocity_qSV_waves_uniform_and_binary_fractal_fractures_300Hz_angle.eps}

\vspace{1cm}

\caption{Phase velocities of qSV waves for a brine saturated sample with 
uniform and binary fractal fractures.}
\label{fig13_veloc_qP_uniform_fractal_fractures}
\end{center}
\end{figure}


\newpage

\begin{figure}[!h]
\begin{center}
\includegraphics[width=16cm]
{fig14_phase_velocity_SH_waves_uniform_and_binary_fractal_fractures_300Hz_angle.eps}

\vspace{1cm}

\caption{Phase velocities of SH waves for a brine saturated sample with 
uniform and binary fractal fractures.}
\label{fig14_veloc_SH_uniform_fractal_fractures}
\end{center}
\end{figure}


\newpage

\begin{figure}[!h]
\begin{center}
\includegraphics[width=16cm]
{fig15_1000_over_qP_uniform_and_binary_fractal_fractures_300Hz_angle.eps}

\vspace{1cm}

\caption{Dissipation factors of qP waves for a brine saturated sample with 
uniform and binary fractal fractures.}
\label{fig15_aten_qP_uniform_fractal_fractures}
\end{center}
\end{figure}


\newpage

\begin{figure}[!h]
\begin{center}
\includegraphics[width=16cm]
{fig16_1000_over_qSV_uniform_and_binary_fractal_fractures_300Hz_angle.eps}

\vspace{1cm}

\caption{Dissipation factors  of qSV waves for a brine saturated sample with 
uniform and binary fractal fractures.}
\label{fig16_aten_qSV_uniform_fractal_fractures}
\end{center}
\end{figure}

\newpage

\begin{figure}[!h]
\begin{center}
\includegraphics[width=16cm]
{fig17_1000_over_q_SH_waves_uniform_and_binary_fractal_fractures_300Hz_angle.eps}

\vspace{1cm}

\caption{Dissipation factors  of SH waves for a brine saturated sample with 
uniform and binary fractal fractures.}
\label{fig17_aten_SH_uniform_fractal_fractures}
\end{center}
\end{figure}

\end{document}

