%\documentclass[doublespacing]{elsart}
\documentclass{elsart}
\usepackage{subeqnarray,amsmath,amsfonts,graphicx}
\usepackage{oldgerm}
\usepackage{rotating}
\bibliographystyle{elsart-num}
\oddsidemargin -.5cm
\journal{Comput. 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}}

\begin{document}

\begin{frontmatter}
\title{\LARGE
Effective viscoelastic medium from  fractured fluid-saturated poroviscoelastic   media. A finite element approach.
}
\author[ba,eu]{Juan E. Santos\corauthref{cor}}
\author[ogs]{Jos\'e M. Carcione}
\address[ba]{CONICET, IGPUBA, Universidad
de Buenos Aires, Argentina}
\address[eu]{Department of Mathematics, Purdue University, 150 N. University  Street, West Lafayette, Indiana, 47907-2067, USA;
santos@math.purdue.edu}
\address[ogs]{OGS, Trieste, Italy}

\corauth[cor]{Corresponding author,  e-mail: santos@math.purdue.edu}

\end{frontmatter}
\newpage
\begin{frontmatter}
\begin{abstract}
An effective viscoelastic medium is derived using harmonic FE numerical experiments in a poroviscoelastic fluid-saturated medium containing a dense set of horizontal fractures. 
\end{abstract}
\begin{keyword}
\sep Poroviscoelasticity   
\sep Finite element methods \sep Effective anisotropic media 

%\PACS 4320.Gp \sep 43.20.Jr
\end{keyword}
\end{frontmatter}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Introduction}
to be written

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{A fractured  Biot's medium}
We consider a porous  solid   saturated by a single phase, 
compressible viscous fluid  and assume 
 that the whole aggregate is isotropic.  
Let $u_{s} = (u_{s,i})$ and 
$\tilde u_{f} = (\tilde u_{f,i}),\, i=1,\cdots,E$ denote the  
averaged displacement vectors  of the solid and fluid phases, 
respectively, where $E$ denotes the Euclidean dimension.
Also let  
\[
u_{f}=\phi(\tilde u_{f}-u_{s}),
\]
be  the average relative fluid  displacement per unit 
volume of bulk material, with $\phi$ denoting the   
effective porosity. Set $u = (u_s,u_f)$ and note that  
\[
\qquad \xi =-\nabla\cdot u_f,
\]
represents the change in fluid content.

Let  $\varepsilon_{ij}(u_s)$  be the   strain 
tensor of the solid.
Also, let $\sigma_{ij},~i,j=1,\cdots,E$, and 
$p_f$   denote  the   stress tensor
of the bulk material  and the fluid pressure, respectively.   
Following \cite{biot62}, the stress-strain relations 
can be written in the form:

\begin{subeqnarray}\label{mod.2}
&&\sigma_{ij}(u)=2 \mu\,\varepsilon_{ij}(u_s)+\delta_{ij}
(\lambda_c \,\nabla\cdot u_s - \alpha ~ K_{av}\, \xi),\slabel{mod.2a}\\
&& p_f(u) =  - \alpha ~ K_{av} \, \nabla\cdot u_s + K_{av} \xi\slabel{mod.2b}.
\end{subeqnarray}
The coefficient $\mu$ is equal to the shear modulus of the bulk 
material, considered 
to be equal to the shear modulus of the dry matrix. Also   
\begin{eqnarray}\label{mod.22}
\lambda_c= K_c - \frac{2}{E} \mu, 
\end{eqnarray} 
with $K_c$ being the bulk modulus of the saturated 
material. Following \cite{santos92} \cite{gassmann51} the 
 coefficients in (\ref{mod.2}) can be obtained from the relations
\begin{eqnarray}\label{mod.3}
\alpha = 1 - \dfrac{K_m}{K_s}, 
\quad  K_{av} = \left(\dfrac{\alpha- \phi}{K_s} 
+ \dfrac{\phi}{K_f}\right)^{-1}\quad K_c = K_m + \alpha^2 K_{av},
\end{eqnarray} 
where $K_s, K_m$ and $ K_f$  denote the bulk modulus of the solid grains 
composing the solid matrix,  the dry matrix and the  the saturant fluid, 
respectively.
The coefficient $\alpha$ is known as the effective stress coefficient of
the bulk material.


\subsection{The equations of motion}
Let $\rho_s$ and $\rho_f$ denote the mass densities of the solid grains 
and the fluid and let 
\begin{eqnarray}\label{density}
&&\rho_b = (1 - \phi) \rho_s  + \phi \rho_f\nonumber
\end{eqnarray}
denote the mass density of the bulk material. 
Let  the positive definite matrix ${\mathcal P}$ and the 
nonnegative matrix ${\mathcal B}$ be defined by  
\begin{eqnarray*}
{\mathcal P} =\begin{pmatrix}
\rho_b I&\rho_f I\\
\rho_f I& m I
\end{pmatrix},
\quad {\mathcal B}=\begin{pmatrix}
0 I&0 I\\
0 I&b I
\end{pmatrix}. 
\end{eqnarray*}
Here $I$ denoted the identity matrix in $R^{E\times E}$. The mass 
coupling coefficient $m$ represents the inertial effects 
associated with
dynamic interactions between the solid and fluid phases, 
while the coefficient 
$b$ includes the viscous coupling effects between  such phases.
They are given by the relations
 
\begin{eqnarray}\label{def_S} 
b= \dfrac{\eta}{k},\qquad m = \dfrac{S \rho_f}{\phi},\qquad 
S = \frac{1}{2} \left(1 + \frac{1}{ \phi}\right),
\end{eqnarray}
where $\eta$ is the fluid viscosity  and $k$ the absolute permeability. $S$ is known
as the structure or tortuosity factor.
Next, let ${\mathcal {L}}(u)$ be 
the second order  differential operator defined by 
\[
{\mathcal {L}}(u) = \left(\nabla\cdot\sigma(u), 
-\nabla p_f(u) \right).
\]
Then if $\o=2\pi f$ is the angular frequency, in the absence of body forces  
 Biot's  equations of motion,     stated in the space-frequency
domain are \cite{biot56b} \cite{biot56c}

\begin{eqnarray}\label{mod10}
-\o^2 {\mathcal {P}} u(x,\o) + i \o 
{\mathcal {B}} u(x,\o) - {\mathcal {L}}(u(x,\o))= 0.
\end{eqnarray}



\subsection{The bounday conditions at a fracture}
Consider a square domain 
$\O=(0,D)^2$  with boundary $\G$ in the $(x_1, x_3)$-plane and 
that we have   one horizontal   
fractures $\G^{(f,1)}$  of length $D$ 
 in our domain $\O$. 
This fracture  divides our domain in two nonoverlapping rectangles $R^{(1)},R^{(2)}$ having as a common side $\G^{(f,1)}$,  so that 
\[
\O = \cup_{l=1}^{2} R^{(l)}.
\] 
Let $\nu_{1,2},\chi_{1,2}$ 
be the  unit outer  normal and a unit tangent (oriented counterclockwise)  
on $\G^{(f,1)}$  from  $R^{(1)}$  to $R^{(2)}$ , 
such that $\{\nu_{1,2}, \chi_{1,2} \}$ is an orthonormal system on $\G^{(f,1)}$.

Let $[u_s], [u_f]$ denote the jumps of the solid and fluid displacement vectors
 at $\G^{(f,1)}$   i.e.,
\[
[u_s] = \left(u_s^{(2)} - u_s^{(1)}\right)|_{\G^{(f,1)}} 
\]
$u_s^{(1)}$ denoting the displacement vector values restricted to $\O^1$, and simularly for $u_s^{(2)}$.

Following Nakawa and Schoemberg, their equation (53) , JASA 2007, we impose 
the following boundary conditions on $\G^{(f,1)}$ (here we use $Z_N$ and $Z_T$ for $\eta_N$ and $\eta_T$ in Nakawa and Schoemberg JASA07)

\begin{eqnarray}
&&[u_s]\cdot \nu_{1,2} = Z_N \sigma(u)\nu_{12} \cdot\nu_{12} 
+  \alpha Z_N p_f(u), \quad \G^{(f,1)}\label{bc_11}\\
&&[u_s]\cdot \chi_{1,2} =  Z_T \sigma(u)\nu_{12} \cdot\chi_{12}\quad \G^{(f,1)}\label{bc_12}\\
&&[u_f]\cdot \nu_{1,2} = - \alpha Z_N \sigma(u)\nu_{12} \cdot\nu_{12} 
- \dfrac{\alpha Z_N}{\widetilde B} p_f(u)\quad \G^{(f,1)}\label{bc_13}.
\end{eqnarray}
Now from \eqref{bc_11} and \eqref{bc_13} we obtain
\begin{eqnarray}
&&\sigma(u)\nu_{12} \cdot\nu_{12} = d_{11} [u_s]\cdot \nu_{1,2} + d_{12} [u_f]\cdot \nu_{1,2}\label{eq_1}\quad \G^{(f,1)}\\
&&- p_f(u) = d_{12} [u_s]\cdot \nu_{1,2} + d_{22} [u_f]\cdot \nu_{1,2}\quad \G^{(f,1)}\label{eq_2},
\end{eqnarray}
where 
\begin{eqnarray}\label{def_fij}
&&f_{11} = \dfrac{\alpha Z_N}{\widetilde B}/\Theta, \quad
 f_{12} = \alpha Z_N/\Theta, \quad f_{22} =  Z_N/\Theta,\\
&&\Theta = \alpha Z_N^2 \left(\dfrac{1}{\widetilde  B} - \alpha\right).\nonumber
\end{eqnarray}
Set 
\begin{eqnarray} \label{bc}
{\mathbf F}=
\begin{pmatrix}
f_{11}&0&f_{12} \\
0&\dfrac{1}{Z_T}&0\\
f_{12}&0&f_{22}
\end{pmatrix} 
\end{eqnarray} 
Then we can write the boundary conditions \label{bc_1}-\label{bc_3} as follows: 
\begin{eqnarray} \label{bc_final}
&&\left(\sigma(u)\nu_{12} \cdot\nu_{12}, \sigma(u)\nu_{12} \cdot\chi_{12}, - p_f(u)\right) \\
&&\hskip1cm = {\mathbf F} \left(
[u_s]\cdot \nu_{1,2}, [u_s]\cdot \chi_{1,2}, [u_f]\cdot \nu_{1,2}\right)^T,\nonumber
\end{eqnarray}
where $^T$ indicates the traspose.

  
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=D\},\nonumber\\
&&\G^B = \{ (x_1, x_3)\in \G: x_3=0\},\quad \G^T = \{ (x_1, x_3)\in \G: x_3=D\}.\nonumber
\end{eqnarray}

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

For obtaining the complex  coefficient $p_{33}$  
let us consider the solution of (\ref{mod10}) with  the boundary condition \eqref{bc_final} together with 
the following boundary conditions 
\begin{eqnarray}\label{bc1a}
&&\sigma(u) \nu \cdot \nu = -\Delta P, \quad (x,y)\in \G^T,\label{bc1b}\\
&&\sigma(u) \nu \cdot \chi = 0, \quad (x,y)\in \G^T,\label{bc1c}\\
&&\sigma(u) \nu \cdot \chi = 0, \quad (x,y)\in \G^L\cup\G^R,\label{bc1d}\\
&&u^s\cdot\nu = 0, \quad (x,y)\in \G^L\cup\G^R,\label{bc1e}\\
&&u^s = 0, \quad (x,y)\in \G^B,\label{bc1f}\\
&&u^f\cdot\nu=0, \quad (x,y)\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{mod10} with the boundary
conditions \eqref{bc1b}-\eqref{bc1g} and \eqref{bc_final}, the vertical displacements  $u^s_y(x,L,\omega)$ on  $\Gamma^T$ allow 
 us to obtain an average vertical displacement $u^{s,T}_y(\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^{s,T}_y(\omega)$, 
which enable us to compute  $p_{33}(\omega)$ by using the relation (\ref{K-dv}). 
The corresponding 
complex compressional velocity is  
\begin{equation}\label{complex_vp}
V_{pc}(\o) = \sqrt{\frac{p_{33}(\o)}{\overline \rho_b}},
\end{equation}
 where $\overline \rho_b$ is the average bulk 
density of the sample.


The following relations allow us to estimate the {\it equivalent}  
compressional phase velocity $V_p(\o)$ and quality factor $Q_p(\o)$ in the form 
\cite{rubino07}:
\begin{eqnarray}\label{vp_qp}  
V_p(\o) = \left[ Re\left(\frac1{V_{pc}(\o)}\right)\right]^{-1},    
\qquad \frac1{Q_p(\o)}=\frac{\text{Im} (V_{pc}(\o)^2)}{\text{Re} (V_{pc}(\o)^2)}
\end{eqnarray}


For obtaining the {\it equivalent} 
complex shear modulus of our fluid-saturated porous medium, 
let us consider the solution of (\ref{mod10}) with \eqref{bc_final} and 
the following boundary conditions 
\begin{eqnarray}\label{s_bc1a}
&&-\sigma(u) \nu  = g, \quad (x,y)\in \G^T\cup\G^L\cup\G^R,\label{s_bc1b}\\
&&u^s = 0, \quad (x,y)\in \G^B,\label{s_bc1f}\\
&&u^f\cdot\nu=0, \quad (x,y)\in \G\label{s_bc1g}, 
\end{eqnarray}
where 
\begin{eqnarray*}
g=\begin{cases} (0, \Delta p),\quad (x,y)\in \G^L,\\
(0, -\Delta p),\quad (x,y)\in \G^R,\\
(-\Delta p,0), \quad (x,y)\in \G^T.\end{cases} 
\end{eqnarray*}
The change in shape of the rock sample allows to recover its {\it equivalent}  
complex shear modulus  $\overline \mu_c(\omega)$ by using the relation
\begin{equation}\label{tg}
tg(\beta\omega))=\frac{\Delta p}{\overline\mu_c(\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_x(x,L,\omega)$ at the top 
boundary $\Gamma^T$ allow  us to obtain, for each frequency,  an average horizontal 
displacement  $u^{s,T}_x(\omega)$ suffered by the 
boundary $\Gamma^T$. This average value allows us to approximate the 
change in shape suffered by
the sample, given by $tg(\beta(\omega))\approx u^{s,T}_x(\omega)/L$, 
which from \eqref{tg}
let us estimate $\overline \mu_c(\omega)$.

The complex shear velocity is given by
\begin{equation}\label{complex_vs}
V_{sc}(\omega) = \sqrt{\frac{\overline\mu_c(\omega)}{\overline{\rho_b}}}
\end{equation}
and the {\it equivalent}  shear phase velocity $V_s(\omega)$ and (inverse) 
quality factor $Q_s(\omega)$ are estimated using 
the relations
\begin{eqnarray}\label{vs_qs}  
V_s(\omega)=
\left[\mbox{Re}\left(\frac1{V_{sc}(\omega)}\right)\right]^{-1},\quad   
\frac{1}{Q_s(\omega)}=\frac{\mbox{Im} (V_{sc}(\omega)^2)}{\mbox{Re} (V_{sc}(\omega)^2)}.
\end{eqnarray}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{A variational formulation}
In order to state a variational formulation for \eqref{mod10} and either 
\eqref{bc1a}-\eqref{bc1g} or \eqref{s_bc1a}-\eqref{s_bc1g} 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)$.
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$.


 Also, let us introduce the spaces 
\[
{\mathcal W}_{11}(\O)= \{v\in [L^2(\O)]^2:  v|_{R^{(l)}}\in [H^1(R^{(l)})]^2 , l=1,2 \,\, v\cdot\nu = 0 \,\,\,  \text{on} \, \,\G^B\cup \G^L\cup \G^R\},
\]

\[
H_0(\opdiv;\O)=\{v\in [L^2(\O)]^2:\, v|_{R^{(l)}}\in H(\text{div}, R^{(l)}) , v\cdot\nu = 0
\,\,\,\text{on} \,\,\G\},
\]
Multiply equation (\ref{mod10}) by $v=\left(v^{s},v^{f}\right)\in
 {\mathcal W}_{11}(\O)\times H_0(\opdiv;\O)$, 
use integration by parts  and apply the boundary conditions \eqref{bc1b},
\eqref{bc1c} and \eqref{bc1d} \eqref{bc_final}
to see that the solution 
of (\ref{mod10}) and ({\ref{bc1a})-\eqref{bc1g} and \eqref{bc_final}  
satisfies {\it the weak form}: 

Set 
\begin{eqnarray}\label{weak_p33}
&&\Lambda(u,v) = -\o^2 \left({\mathcal P} u,v\right) + i \o \left({\mathcal B} u,v\right)\\ 
&&+ \sum_{l=1}^2\sum_{s,t}\left(\tau_{st}(u),\e_{st}(v^{s})\right)_{R^{(l)}}
 - \left(p_f(u),\nabla\cdot v^{f})\right)_{R^{(l)}}\nonumber\\
%&& \left\langle
%{\mathbf F} 
&&\left\langle{\mathbf F}\left(
[u_s]\cdot \nu_{1,2}, [u_s]\cdot \chi_{1,2}, [u_f]\cdot \nu_{1,2}\right)^T,
\left(
[v_s]\cdot \nu_{1,2}, [v_s]\cdot \chi_{1,2}, [v_f]\cdot \nu_{1,2}\right)\right\rangle\nonumber\\
&&= \left\langle \Delta P, v^{s}\cdot\nu\right\rangle_{\G^T}, \qquad  
\forall v=\left(v^{s},v^{f}\right)\in {\mathcal W}_{11}(\O)\times H_0(\opdiv;\O)\nonumber
\end{eqnarray}

%&&\qquad =  -\o^2 \left({\mathcal P} u,v\right) + i \o \left(b u^f,v^f\right) 
%+  \sum_{l=1}^{2}\left({\mathbf E} ~ \tilde \epsilon(u),\tilde %\epsilon(v)\right)_{R^{(l)}}\nonumber\\
% \nonumber
%\end{eqnarray}
Note  that  in \eqref{weak_p33},  we can write 
\begin{eqnarray}
&&\sum_{l=1}^2\sum_{s,t}\left(\tau_{st}(u),\e_{st}(v^{s})\right)_{R^{(l)}}
 - \left(p_f(u),\nabla\cdot v^{f})\right)_{R^{(l)}}\label{strain}
%\\
%&& 
= \sum_{l=1}^{2}\left({\mathbf E} ~ \tilde \epsilon(u),\tilde \epsilon(v)\right)_{R^{(l)}}
\end{eqnarray}
In \eqref{strain}, the   complex 
matrix ${\mathbf E(\o)} = {\mathbf E(\o)}_R + i {\mathbf E(\o)}_I$ and the column vector  $\tilde\varepsilon((u))$
are defined  by
\begin{eqnarray*}
{\mathbf E}=\begin{pmatrix}
\lambda_c+2\mu&\lambda_c&\alpha K_{av}&0 \\
\lambda_c&\lambda_c+2\mu&\alpha K_{av}&0\\
\alpha K_{av}&\alpha K_{av}&K_{av}&0\\
0&0&0&4 \mu
\end{pmatrix}&,&\,\,
\tilde\varepsilon(u) = \begin{pmatrix}
\e_{11}(u^s) \\ \e_{22}(u^s) \\ \nabla\cdot u^f\\ \e_{12}(u^s)
\end{pmatrix}.
\end{eqnarray*}
It will be assumed that the real part ${\mathbf E}_R(\o)$ is 
positive definite since in the
elastic limit it is associated with the strain energy density. 
Furthermore, the imaginary parts  ${\mathbf E}_I(\o)$ are  assumed to be positive definite because of
the restriction imposed on our system  by the first and second laws of 
thermodynamics. 

\vskip1cm

{\bf Jose:  aqui hay que asumir que ${\mathbf E}_I(\o)$ es definida positiva, o sea poroviscoelasticidad, para poder demostrar que la solucion de los $p_{IJ}$
es unica. Despues en los exprimentos, se hara como en el paper en GJI.

Lo que importa es que creas el argumento que lleva a las boundary conditions, y si tenes ganas, escribi la demostracion de  que ese termino $\Theta$ es positivo en el caso general.  Lo interesante es que   la weak form queda muy parecida a la del caso viscoelastico fracturado. 

demas, hay que generalizar las B. C. para que la matriz  $F$ tenga parte imaginaria (positiva), asi tenemos atenuacion en las fracturas, como en el sfrac.pdf. para el $\eta_T$ es facil, es lo mismo del sfrac.pdf. 

todavia no escribi la matematica de ese paper nuesto en GJI, quedo alli en la cola. 
La de este contiene  a esa, asi que ahora a escribir los teoremas antes que me agarre la licuadora de los K.....  } 

\begin{thebibliography}{}
%\baselineskip=2\baselineskip

\bibitem{biot56b}
M. A. Biot,``Theory of propagation of elastic waves in a fluid-saturated 
porous solid. I. Low frequency range'', {J. Acoust. Soc. Am.},{\bf 28},
168--171 (1956).


\bibitem{biot56c}
M. A. Biot,``Theory of propagation of elastic waves in a fluid-saturated 
porous solid. II. High frequency range'', {J. Acoust. Soc. Am.},{\bf 28},
179--191 (1956).

\bibitem{biot62} M. A. Biot, ``Mechanics of deformation and acoustic
propagation in porous media,'' {J. Appl. Phys.} {\bf 33}, 1482--1498  
(1962).

\bibitem{cadoret_95} T. Cadoret and  D. Marion and B. Zinszner, 
``Influence of frequency and fluid distribution on elastic waves velocities in
partially saturated limestones'',{J. Geophysical Research}, {\bf 100},
9789--9803, 1995.



\bibitem{car-pic06}
J. M. Carcione and  S. Picotti,
``P-wave seismic attenuation by slow-wave diffusion: Effects of 
inhomogeneous rock properties  '', 
{Geophysics} {\bf 71}, 01-08 (2006).

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

\bibitem{jdsss}
J. Douglas Jr. and  J. E. Santos and  D. Sheen and L. Bennethum,
``Frequency domain treatment of one--dimensional scalar waves'',
{Math. Models Methods Appl. Sci. }, {\bf 3}, 171-194 (1993).

\bibitem{dfp_97}
J. Douglas~Jr. and F. Furtado and F. Pereira,
``The statistical behavior of instabilities in immiscible displacement subject
to fractal geology'',
{Mathematical Modeling of Flow through Porous Media, eds. A. Burgeaut, C. 
Carasso, S. Luckhaus and A. Mikelic, World Scientific, Singapore}, 115--137, 
 (1997).



\bibitem{dss01}
J. Douglas~Jr. and J. E. Santos and D. Sheen,
``Nonconforming Galerkin methods for the Helmholtz equation'',
{Numer. Methods for Partial Diff. Equations}, {\bf  17}, 475-494 (2001).


\bibitem{duttayode}
N. C. Dutta, H. Od\'e, ``Attenuation and dispersion of compressional waves in
fluid-filled porous rocks with partial gas saturation (White model). Part I: Biot
theory''
 {Geophysics}, {\bf 44 (11)},
1777-1788 (1979).

\bibitem{fc} A. Frankel and 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.} {\bf  91} 6465 - 6489,  (1986).

\bibitem{gassmann51}
F. Gassmann, ``\"Uber die elastizit\"at por\"oser medien'' (``On the elasticity   of porous media''), { Vierteljahrsschrift 
der Naturforschenden Gessellschaft in Zurich} {\bf 96}, 1--23 (1951).

\bibitem{luis_mathgeo_04}
L. Guarracino and J. E. Santos, 
Stochastic modeling of variably saturated flow in fractal porous media , 
 Mathematical Geology {\bf 36 (2)} 239--260, (2004) 

\bibitem{hpc03}
H.B. Helle, N. H. Pham and J.M. Carcione, 
Velocity and attenuation in partially saturated rocks - 
Poroelastic numerical experiments, 
Geophysical Prospecting {\bf 51} 551-566 (2003).

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


\bibitem{mp07}
Y. Mason  and  S. Pride, 
``Poroelastic finite difference modeling of seismic attenuation and dispersion due to mesoscopic-scale heterogeneity'', 
Journal of Geophysical Research, {\bf 112}, B03204. 


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

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

\bibitem{ptb02}
S. Pride  and E. Tromeur and J. G. Berryman, 
``Biot slow-wave effects in stratified rock''
{Geophysics}, {\bf 67}, 271-281, (2002).



\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).

\bibitem{rubino07}
J. G. Rubino, J. E. Santos, S. Picotti and J. M. Carcione,
`Simulation of upscaling effects  due to wave-induced fluid flow in Biot media using the finite element method��, {Journal of Applied Geophysics}, {\bf 62}, 193-203, (2007).

\bibitem{santos92}
J. E. Santos, J. M. Corbero, C. L. Ravazzoli, J. L. Hensley, ``Reflection and transmission
coefficients in fluid-saturated porous media,'' {J. Acoust. Soc. Am.}{\bf 91 (4)},
 1911-1923 (1992).
 
\bibitem{santos_comp_geo_05} J. E. Santos. C. L. Ravazzoli, P. M. Gauzellino and J.
M. Carcione,
Numerical simulation of ultrasonic waves in reservoir rocks with 
 patchy saturation 
and fractal petrophysical properties, {Computational Geosciences} {\bf 9} (2005) 1-27.



\bibitem{white75}
J. E. White, ``Computed seismic speeds and attenuation in 
rocks with partial gas saturation'', {Geophysics}, {\bf 40}, 224-232, (1975).

\bibitem{white75b}
J. E. White and N. G. Mikhaylova and F. M. Lyakhovitskiy,
``Low-frequency seismic waves in fluid-saturated layered rocks'',
{Izvestija Academy of Siences USSR, Physics of Solid Earth}, {\bf 10}, 654-659, (1975).




\end{thebibliography}

\newpage
\noindent{FIGURE CAPTIONS}
\vskip1cm
{\bf Figure 1:} P-wave phase velocity obtained from the compressibility test (dots) 
and using White's theory (line) for frequencies lying between $0$  and $100$ Hz.

\vskip .5cm
{\bf 
Figure 2:} P-wave inverse quality factor obtained from the compressibility test 
(dots) and using White's theory (line) for frequencies lying between $0$  and $100$ Hz.

%**************** 
\vskip.5cm
{\bf 
Figure 3:} Gas-water distribution for a given realization. Black zones
correspond to pure gas saturation and the white ones to pure
water saturation. The overall gas saturation is $0.1$.

\vskip.5cm
{\bf 
Figure 4:} Normalized fluid-pressure amplitude for the fluid distribution
shown in Figure 3. The excitation frequency is $190$ Hz.

\vskip.5cm
{\bf 
Figure 5:} Averaged variance of the compressional phase velocity as  function of
the total number of realizations.

\vskip.5cm
{\bf 
Figure 6:} Effective compressional phase velocity as  function of
frequency (solid lines). Dotted lines indicate the corresponding
standard deviations.

\vskip.5cm
{\bf 
Figure 7:} Effective compressional inverse quality factor as  function
of frequency.  Dotted lines indicate the corresponding
standard deviations.


%****************** 

\vskip.5cm
{\bf 
Figure 8:} Distribution of shale and sandstone. Black zones correspond to
pure shale, and the white ones to pure sandstone.

\vskip.5cm
{\bf 
Figure 5:} Equivalent P-wave phase velocity as  function of frequency.

\vskip.5cm
{\bf 
Figure 10:} Equivalent compressional inverse quality factor as  function
of frequency.

\vskip.5cm
{\bf 
Figure 11:} Equivalent shear  inverse quality factor as  function
of frequency.


\end{document}
\begin{figure}
%\begin{center}
%\includegraphics[scale=.3]{fig1a.eps}
%\includegraphics[scale=.4]{fig1b.eps}
\caption{\label{fig1} The figure on the left shows the degrees of freedom of a NC-stencil. 
However, in order to perform a dispersion analysis of the NCFEM method, a stencil without
preferred direcions is needed. The figure on the right shows a scheme of how it is 
obtained. In both figures the filled circles represent the degrees of freedom corresponding
to the solid,  the circle to the $x_1$-component of the fluid and crosses to the $x_2$-component of the fluid. }
%\end{center}
\end{figure}
%
\newpage
\pagestyle{empty}
\begin{figure}
%\unitlength 1 mm
%\begin{picture}(150,50)(0,0)
%\put(-10,0){\resizebox{.5\textwidth}{!}{\includegraphics{fig2a.eps}}}
%\put(65,0){\resizebox{.5\textwidth}{!}{\includegraphics{fig2b.eps}}}
%\put(140,0){\resizebox{.5\textwidth}{!}{\includegraphics{fig2c.eps}}}
%\end{picture}
\caption{\label{fig2} Non-viscous case. Dimensionless phase velocities 
for the fast (left), slow (middle) and shear waves (right) as a function 
of the reciprocal of the 
number of points per wavelegth for different propagation directions.}
\end{figure}
%
\begin{figure}
%\unitlength 1 mm
%\begin{picture}(150,50)(0,0)
%\put(-15,0){\resizebox{.4\textwidth}{!}{\includegraphics{fig3a.eps}}}
%\put(70,0){\resizebox{.4\textwidth}{!}{\includegraphics{fig3b.eps}}}
%\put(155,0){\resizebox{.4\textwidth}{!}{\includegraphics{fig3c.eps}}}
%\end{picture}
\caption{\label{fig3} Viscous case. Dimensionless phase velocities for the 
fast (left), slow (middle) and shear waves (right), as a
function of the frequency, for different numbers of points per wavelength. 
Only one propagating direction is depicted; for most cases this
corresponds to the worst one.}
\end{figure}
%
\begin{figure}
%\unitlength 1 mm
%\begin{picture}(150,50)(0,0)
%\put(-15,0){\resizebox{.4\textwidth}{!}{\includegraphics{fig4a.eps}}}
%\put(70,0){\resizebox{.4\textwidth}{!}{\includegraphics{fig4b.eps}}}
%\put(155,0){\resizebox{.4\textwidth}{!}{\includegraphics{fig4c.eps}}}
%\end{picture}
\caption{\label{fig4} Same as Fig.\ref{fig3} for the group velocity case.}
\end{figure}
%
\begin{figure}
%\unitlength 1 mm
%\begin{picture}(150,50)(0,0)
%\put(-15,0){\resizebox{.4\textwidth}{!}{\includegraphics{fig5a.eps}}}
%\put(70,0){\resizebox{.4\textwidth}{!}{\includegraphics{fig5b.eps}}}
%\put(155,0){\resizebox{.4\textwidth}{!}{\includegraphics{fig5c.eps}}}
%\end{picture}
\caption{\label{fig5} Same as Fig.\ref{fig3} for the attenuation.}
\end{figure}
%
\begin{figure}
%\unitlength 1 mm
%\begin{picture}(150,50)(0,0)
%\put(-15,0){\resizebox{.4\textwidth}{!}{\includegraphics{fig6a.eps}}}
%\put(70,0){\resizebox{.4\textwidth}{!}{\includegraphics{fig6b.eps}}}
%\put(155,0){\resizebox{.4\textwidth}{!}{\includegraphics{fig6c.eps}}}
%\end{picture}
\caption{\label{fig6} Viscous case. 
Dimensionless phase velocities for the fast (left), slow (middle) and 
shear (right) waves, as a function of the reciprocal of the number of points per
wavelength at a fixed frequency for different propagation directions.}
\end{figure}
%
\begin{figure}
%\unitlength 1 mm
%\begin{picture}(150,50)(0,0)
%\put(-15,0){\resizebox{.4\textwidth}{!}{\includegraphics{fig7a.eps}}}
%\put(70,0){\resizebox{.4\textwidth}{!}{\includegraphics{fig7b.eps}}}
%\put(155,0){\resizebox{.4\textwidth}{!}{\includegraphics{fig7c.eps}}}
%\end{picture}
\caption{\label{fig7} Viscous case. Dimensionless group velocities for the fast (left), 
slow (middle) and shear (right) waves, as a function of the reciprocal of the number of points per
wavelength at a fixed frequency for different propagation directions.}
\end{figure}
%
\begin{figure}
%\unitlength 1 mm
%\begin{picture}(150,50)(0,0)
%\put(-15,0){\resizebox{.4\textwidth}{!}{\includegraphics{fig8a.eps}}}
%\put(70,0){\resizebox{.4\textwidth}{!}{\includegraphics{fig8b.eps}}}
%\put(155,0){\resizebox{.4\textwidth}{!}{\includegraphics{fig8c.eps}}}
%\end{picture}
\caption{\label{fig8} Viscous case. Dimensionless attenuations for the 
fast (left), slow (middle) and shear (right) waves, as a function of the reciprocal of the number of points per
wavelength at a fixed frequency for different propagation directions.}
\end{figure}
%\nonumsection{REFERENCES }
  





















