%\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}
\newcommand{\thmref}[1]{Theorem~\ref{#1}}
\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}}
%       Theorem environments

%\newenvironment{proof}{\noindent{\bf Proof:\ }}{\hspace*{\fill} \QED} 

\newtheorem{theorem}{Theorem}
\newenvironment{proof}{\noindent{\it Proof:\ }}{\hspace*{\fill} {\it QED}}



%\newenvironment{endproof}{\hspace*{\fill} {\bf QED}}
%[section]

\begin{document}

\begin{frontmatter}
\title{\LARGE
Finite Element Procedures for 2D and 3D Seismoelectric Modeling
}
\author[obs,eu]{Juan E. Santos\corauthref{cor}}
\author[obs]{Fabio I. Zyserman}
\author[obs]{Patricia M. Gauzellino}
\address[obs]{CONICET, Fac. de Cs. Astron\'omicas y Geof\'{\i}sicas, Universidad
Nacional de La Plata,
 Paseo del Bosque s/n, B1900FWA - La Plata, Argentina}
\address[eu]{Department of Mathematics, Purdue University, 150 N. University  Street, West Lafayette, Indiana, 47907-2067, USA;
santos@math.purdue.edu}


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

\end{frontmatter}
\newpage
\begin{frontmatter}
\begin{abstract}


This work  presents a collection continuous and discrete time  finite element 
procedures for  2D and 3D seismoelectric modeling in poroelastic
fluid-saturated media.  The  model  involves  the simultaneous  
 solution  of Biot's equations of motion and Maxwell 
equations, coupled via an electrokinetc coefficient, employing  
 absorbing boundary conditions at the artificial boundaries.
The case of compressional and vertically polarized shear waves 
 coupled with the transverse magnetic polarization  (PSVTM-mode)    
is analyzed in detail, including the  derivation of {\it a priori} error 
estimates on the 
 continuous-time finite element procedure and stability results for the two discrete
 time procedures formulated, both for  rectangular and triangular elements. Maxwell
 equations are discretized in space  using the mixed finite element spaces of Nedelec, while 
 while for Biot' s equations  we used a nonconforming element for 
 each component of the solid displacement and the vector part of the
 Raviart-Thomas-Nedelec of zero order for the fluid  displacement vector.   
 Later the corresponding  finite element spaces and continuous and discrete-time 
 procedures for the full 3D case are formulated and analyzed both for simplicial and  parallelepiped
 elements.
\end{abstract} 
\begin{keyword}

Seismoelectric Modeling\sep Poroelasticity \sep Electromagnetics  
\sep Finite element methods \sep  

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Introduction}

\vskip1cm
{\bf cambiar es identico al paper en NMPDE}
\vskip1cm

Seismic waves propagating through
near-surface layers of the Earth  may  induce 
electromagnetic disturbances 
that can be measured at the surface (seismoelectric effect)\cite{pride96,mikha97,mikha2000}. Also,  recent tests
suggest that the reciprocal process, {\it i.e.}  surface measurable  
acoustic disturbances induced by 
 electromagnetic fields (electroseismic effect),  is also possible
\cite{thompson05a,hornbostel07}.\\
In order to explain these phenomena, Thompson and  Gist \cite{thompson93}
and Pride \cite{pride_94} suggested that they are generated by an
electrokinetic coupling
mechanism which  can be
shortly explained as follows \cite{block06,haines_pride_06}.  Within a fluid saturated porous medium there
exists a nanometer-scale separation of electric charge in which a bound
charge existing on the surface of the solid matrix (normally of negative sign)
is balanced by adsorbed positive ions of the surrounding fluid, setting an immobile layer.
Further from the surface there exists a distribution of mobile counter ions,
forming the so called diffuse layer. The effective thickness of this double
layer is of about 10 nm. When an electric field is applied to
this system, the ions in the diffuse layer move, dragging the pore fluid
along with it because of the viscous traction. This is known as electro-osmosis
and is responsible for the electroseismic phenomena.
On the other hand, the reciprocal situation arises when an applied pressure gradient creates
fluid flow and hence, an ionic convection current, which in turn produces
an electric field. This is known as electrofiltration and is responsible for the
so-called seismoelectric phenomena.\\
Using a volume averaging approach, Pride \cite{pride_94} derived a set of
equations   describing  both electroseismic and seismoelectric effects in
electrolyte-saturated porous media. In these equations the coupling mechanism
acts through the (generally frequency dependent) electrokinetic coupling
coefficient $L(\o)$. When this coefficient is set to zero, Pride's
set of equations turns to the uncoupled Maxwell's and Biot's equations,
describing the latter mechanical wave propagation in a fluid
saturated porous medium \cite{biot56b,biot56c}.\\
There exist already some works implementing different numerical methods
to solve the set of equations modeling both mentioned processes.
Haartsen and Pride \cite{haartsen_pride_97}, Han and Wang \cite{han_01}, Pain et. al \cite{pain2005}, Haines and Pride
\cite{haines_pride_06}  and White \cite{bwhite_05} and White and Zhou \cite{bwhite_06} have
proposed several different approaches to numerically study these phenomena.\\


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%   SECTION    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{The  differential model}
Consider  a 3D-rectangular  domain $\O = \O_a \cup \O_p$ 
 where $\O_a$ and $\O_p$ 
are  associated with the air and subsurface poroelastic  (disjoint)  parts 
of $\O$, respectively. We will assume that (in cartesian 
coordinates $(x_1, x_2, x_3)$) all physical quantities 
  describing our domains  $\O_a$ and $\O_p$ 
are independent of the $x_2$-direction (i.e., $x_2$ is the symmetry axis) 
and consider  a  seismic line source  in the $x_2$-direction.  Including the air region
allows us to include in implicit fashion the boundary conditions for the
electromagnetic fields  at the air-poroelastic interface.

Let us denote by $E, H$ 
to the electric and magnetic fields in $\O$, respectively,  
and by $u^s, u^f$ to the solid and relative fluid 
displacement vectors in $\O_p$. 
Under the above  symmetry asumption, this source term induces 
 electric and magnetic  fields of the form  
$(E_1(x_1, x_3, t),0,E_3(x_1, x_3, t))$,  
$(0,H_2(x_1, x_3, t),0)$, respectively,  and  solid and relative fluid displacements 
of the form $u^s = (u^s_1(x_1, x_3, t), 0 ,
u^s_3(x_1, x_3, t))$ and  $u^f = (u^f_1(x_1, x_3, t), 0 , u^f_3(x_1, x_3, t))$,
respectively. Consequently only 
 compressional and vertically polarized shear  seismic waves (PSV-waves) are 
 generated. This is a 2D model known as a PSVTM-mode.
 
Let us identify the 3D vectors $(E_1(x_1, x_3, t), 0, E_3(x_1,x_3,t))$ and 
$(0,H_2(x_1, x_3, t)),0)$  with the 2D vector 
$(E(x_1, x_3, t) = (E_1(x_1, x_3, t), E_3(x_1, x_3, t))$ and the scalar 
$H_2(x_1, x_3, t)$, respectively. 
Then  recall that 
\begin{eqnarray*}
\text{curl} ~ H_2 = \left(-\frac{\partial H_2}{\partial x_3}, 
\frac{\partial H_2}{\partial x_1}\right),\quad 
\text{curl}~ E =  \frac{\partial E_1}{\partial x_3} - \frac{\partial E_3}{\partial x_1}.
\end{eqnarray*} 

Also, let us identify our 3D-rectangular  domain $\O$ with   the  
 2D-rectangular domain $\O\cap \{y=0\}$, so that $\O$ 
 is the union of the disjoint rectangular subdomains $\O_a$ and $\O_p$. 
 Let $\G$ denote the boundary of $\O$ and let 
 $\G_{a,p} = \overline \O_a \cap \overline \O_p$ denote the free surface. 
Also let $\G_a = \p\O_a \setminus \G_{a,p}$, $\G_p = \p\O_p \setminus \G_{a,p}$ 
denote the artificial boundaries of $\O_a$ and $\O_p$, respectively. 

Following \cite{pride_94,haines_pride_06}, for 2D seismoelectric modeling 
the electric and  magnetic fields $E$ and $H$ 
and the displacement vectors $u^s$  and 
$ u^f$ satisfy the 
 coupled electromagnetic-poroelastic  equations, 
stated in the space-time domain as follows:

\begin{eqnarray} 
&& \varepsilon \frac{\p E}{\p t}  +  \sigma E - \text{curl}~ H_2 + L_0 \frac{\eta}{\kappa_0} \frac{\p
u^f}{\p t}  = 0,\quad \O,\label{mod.1b}\\
&&\text{curl} ~ E +   \frac{\p  H_2}{\p t} = 0, \quad \O, 
\label{mod.1a}\\
&&  \rho_b \frac{\p^2 u^s}{\p t^2}  + \rho_f \frac{\p^2 u^f}{\p t^2}  
 -  \nabla\cdot \tau(u) = F^{(s)},\quad \O_p,
\label{mod.1c}\\
&& \rho_f \frac{\p^2 u^s}{\p t^2} +  m  \frac{\p^2 u^f}{\p t^2}  + 
  \frac{\eta}{\kappa_0} \ \frac{\p u^f}{\p t} u^f  
 + \nabla p_f = F^{(f)},\quad \O_p,\label{mod.1da}\\
&&\tau_{lm}(u)=2 G\,\varepsilon_{lm}(u^s)+\delta_{lm}
\left(\lambda_c \,\nabla\cdot u^s + \alpha K_{av}\, \nabla\cdot u^f\right),\quad
\O_p,\label{mod.1e}\\
&& p_f(u) =  - \alpha K_{av} \, \nabla\cdot u^s - K_{av}  \nabla\cdot u^f, \quad \O_p.\label{mod.1f}
\end{eqnarray}

In the equations above $u = (u^s, u^f)$ and $\tau_{lm}(u)$ 
is the stress tensor of the bulk material
and $p_f(u)$ the fluid pressure, while $\varepsilon_{lm}(u^s)$  
denotes the strain tensor of the solid frame. The coefficients in the stress strain
relation are determined as follows. The coefficient $G$  
is equal to the elastic shear modulus of the 
 dry matrix.  Also, 
\begin{equation}\label{deflambdac_aster}
\lambda_c= K_c - G, 
\end{equation}
with $K_c$ being the bulk modulus of the saturated 
material. The 
 coefficients in \eqref{mod.1e}-\eqref{mod.1f}  can be obtained from the relations \cite{gassmann51}, \cite{santos92}
\begin{eqnarray}\label{mod.3}
&&\alpha = 1 - \dfrac{K_m}{K_s}, 
\qquad  K_{av} = \left[\dfrac{\alpha- \phi}{K_s} + \dfrac{\phi}{K_f}\right]^{-1}\\
&&K_c = K_m + \alpha^2 K_{av},\nonumber
\end{eqnarray} 
where $\phi$ denotes the effective porosity and 
$K_s, K_m$ and $ K_f$  denote the bulk modulus of the solid grains 
composing the solid matrix,  the dry matrix and the  saturant fluid, 
respectively.  


Also, $\varepsilon$ is the electric permitivity,  $\mu$  the magnetic permeability  
and $\sigma$ the conductivity, while 
 $F^{(s)}$, $F^{(f)}$ are the external 
 seismic sources.  Furthermore,   
\begin{equation}\label{robulk} 
\rho_b = \phi \rho_f + (1-\phi)\rho_s, 
\end{equation} 
where  $\rho_s$ and  $\rho_f$ denote the mass densities of the solid 
grains composing the solid matrix and the saturant fluid. On the other hand,  
$\eta$  is the fluid viscosity,  $\kappa_0$ the   permeability and $m$ 
is the mass coupling
coefficient between the solid and fluid phases in $\O_p$. 
Following \cite{berry}, $m$  can be written in the form

\begin{equation}\label{def_m} 
m = \frac{\alpha_{\infty} \rho_f}{\phi}, 
\end{equation}  
with $\alpha_{\infty} $ being the formation tortuosity.
 
The positive coupling coefficient   
$L_0$  is defined as  \cite{haartsen_pride_97} by
\begin{eqnarray}\label{def_L0}
L_0 = - \frac{\phi}{\alpha_{\infty}} \frac{\epsilon_0 k_f \zeta}{\eta} \left( 1 - 2 \alpha_{\infty} \frac{\widetilde
d}{\Lambda}\right), 
\end{eqnarray}
with  $\zeta = 0.008 + 0.026 \text{log}_{10}(C_e)$ denoting 
 the zeta potential and $C_e$ being the
electrolyte  molarity. In \eqref{def_L0}  $\epsilon_0$ and $k_f$ are   the vacuum and   fluid permitivities and 
\begin{equation}\label{def_debye}
\widetilde d = \frac{\epsilon_0 k_f k_B T}{e^2 {\it z}^2 N_{ic}} 
\end{equation}
is the Debye length in meters. in \eqref{def_debye} $e$ is 
the electronic charge,  $k_B$ is the Boltzman constant, T is 
the absolute temperature (so that $k_B T$ is the thermal energy) {\it z} is 
the ionic valence and $N_{ic}$ the ionic concentration in ions per meters cubed.
 
To solve equations \eqref{mod.1b}-\eqref{mod.1f} in our 2D domain $\O$  we need a
collection of boundary conditions.  
Let $\G$ denote the boundary of $\O$ and let 
 $\G_{a,p} = \overline \O_a \cap \overline \O_p$ denote the free surface. 
Also let $\G_a = \p\O_a \setminus \G_{a,p}$, $\G_p = \p\O_p \setminus \G_{a,p}$ 
denote the artificial boundaries of $\O_a$ and $\O_p$, respectively. 
  Also, if  $\G_s$ is either an inner interface in $\O$ or a part of the
  boundaries $\G, \G_p$ or $\G_{a,p}$, set 
\begin{subeqnarray}\label{mod8}
{\mathcal G}_{\G_s}(u) &=& \bigg(\tau(u) \nu \cdot \nu, \tau(u)\nu\cdot \chi,p_f(u)\bigg)^t,\slabel{mod8a}\\
S_{\G_s}(u) &=& \left(u^s\cdot\nu,u^s\cdot\chi,
 u^f\cdot\nu\right)^t,\slabel{mod8b}
\end{subeqnarray} 
where , where $^t$ denotes the transpose,  $\nu$ is  the unit outer normal 
on $\G_s$  and 
$\chi$ is a unit tangent on $\G_s$ oriented counterclockwise.  




Then, for $\o>0$ consider the solution of  \eqref{mod.1b}-\eqref{mod.1f}
with  the  absorbing boundary conditions  \cite{dsheenc}, \cite{santos_imajna}
\begin{eqnarray}
&& - \varepsilon^{1/2}  E\cdot \chi + H_2 = 0, \quad \text{on} \quad  \G,\label{modf.5}\\
&& -{\mathcal G}_{\G_p}(u)  =  {\mathcal D}S_{\G_p}(\frac{\p u}{\p t}), 
\quad \text{on}\quad \G_p,\label{modf.6}
\end{eqnarray}
the free surface condition
\begin{eqnarray}
&& -{\mathcal G}_{\G_p}(u)  = 0, 
\quad \text{on}\quad \G_{a,p},\label{modf.6a}
\end{eqnarray}
and the initial conditions
\begin{eqnarray}
&&E(x_1,x_3,t=0) = E_0, \quad H_2(x_1,x_3,t=0) = H_{2,0} 
\label{icond}\\
&& u^s(x_1,x_3,t=0) = u^s_0, \quad u^f(x_1,x_3,t=0) = u^f_0,\nonumber\\
&& \frac{\p u^s}{\p t}(x_1,x_3,t=0) = u^s_1, \quad \frac{\p u^f}{\p
t}(x_1,x_3,t=0) = u^f_1.
\nonumber
\end{eqnarray}

The  matrix ${\mathcal D}$ in \eqref{modf.6} is defined as:
${\mathcal D} 
= {\mathcal R}^{\frac12} {\mathcal S}^{\frac12}{\mathcal R}^{\frac12}
$, where ${\mathcal S } 
= {\mathcal R}^{-\frac12} {\mathcal  M}^{\frac12}{\mathcal R}^{-\frac12}
$
and 
\begin{eqnarray}
{\mathcal R} =\begin{pmatrix}
\rho_b &0&\rho_f \\
0&b&0\\
\rho_f&0&\zeta 
\end{pmatrix},
\quad {\mathcal M}=\begin{pmatrix}
\lambda_c + 2 G &0 &\alpha ~ K_{av}\label{def_RM}\\
0 &G&0\\
\alpha ~ K_{av}&0 & K_{av}
\end{pmatrix},
\end{eqnarray}
where 
\[
b=\rho_b - \frac{(\rho_f)^2}{m}.  
\]
{\emph Remark}: Note that since $\alpha_{\infty}\ge 1$, the matrix ${\mathcal R}$ 
is positive definite. Also, we will requiere that the following conditions 
 be satisfied by the coefficients defining the matrix ${\mathcal M}$:
\begin{subeqnarray}
G>0\slabel{cond_1},\\
\lambda_c + G - \alpha^2 K_{av} >0\slabel{cond_2},\\
K_{av} >0\slabel{cond_3}.
\end{subeqnarray}
Conditions \eqref{cond_1}, \eqref{cond_2} and \eqref{cond_3} are necessary 
and sufficient conditions for the matrix ${\mathcal  M}$ to be positive definite. 
 In particular, the condition \eqref{cond_2}  imposes that  the  inverse of the  
jacketed compressibility coefficient be strictly positive, see \cite{biot62}.  
As a consequence of the positive definitess of the matrices ${\mathcal R}$
 and ${\mathcal M}$, the matrix ${\mathcal D}$ is also positive definite.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%   SECTION    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
\section{A weak formulation}
For  $X\subset \mathbb R^d, d=1,2,3$ with boundary $\partial X$, 
let $(\cdot,\cdot)_X$ denote the complex $L^2(X)$ inner product for scalar, 
vector, or matrix valued functions.
Also, for $s\in\mathbb R$,
$\|\cdot\|_{s,X}$  will denote  
the usual norm  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$. Set 
\begin{eqnarray*}
&&H(\text{curl},\O)=\{\psi\in(L^2(\O))^2: \text{curl} \psi\in L^2(\O)\},\\
&&H(\text{div},\O_p)=\{\psi\in(L^2(\O_p))^2: \nabla\cdot \psi\in L^2(\O_p)\},\\ 
&&H^1(\text{div},\O_p)=\{\psi\in(H^1(\O_p))^2: \nabla\cdot \psi\in H^1(\O_p)\},
\end{eqnarray*}
provided with the natural norms
\begin{eqnarray*}
&&\|\psi\|_{H(\text{curl},\O)}=(\|\psi\|^2_0+\| \text{curl} \psi\|^2_0)^{\frac12},\\
&& \|\psi\|_{H(\text{div},\O_p)}=(\|\psi\|^2_{0,\O_p}+\| \nabla\cdot
\psi\|^2_{0,\O_p})^{\frac12},\\ 
&&\|\psi\|_{H^1(\text{div},\O_p)}=(\|\psi\|^2_{1,\O_p}+\| \nabla\cdot
\psi\|^2_{1,\O_p})^{\frac12}.
\end{eqnarray*}
Recall the integration by parts formulas \cite{gira,dsheenb}
\begin{eqnarray}\label{by_parts1}
&&(\psi, \text{curl} \varphi) - (\text{curl}~ \psi,\varphi)=\la\psi\cdot\chi,\varphi\ra,\,\quad
\forall\psi\in H(\text{\text{curl}},\O),\quad \varphi\in H^1(\O),\\
&&(\nabla\cdot \psi, \varphi) 
+ (\psi,\nabla \varphi)=\la\psi\cdot\nu,\varphi\ra,\,\quad
\forall\psi\in H(\text{div},\O),\quad \varphi\in H^1(\O).\label{by_parts2}
\end{eqnarray}
To obtain a variational formulation, 
test \eqref{mod.1a} against $\varphi\in L^2(\O)$  and test \eqref{mod.1b} 
against $\psi\in H(\text{curl},\O)$ and use 
the integration by parts formula \eqref{by_parts1} and the 
boundary condition \eqref{modf.5} to obtain the equations
\begin{eqnarray}
&& (\varepsilon \frac{ \p E}{\p t}, \psi)
 + (\sigma E, \psi) -(H_2, \text{curl}~ \psi)  
 + \left( L_0 \frac{\eta}{\kappa_0} \frac{\p u^f}{\p t}, \psi\right)_{\O_p}
 \label{eq_13}\\
&&\hskip3cm + \left\langle \left( \frac{\varepsilon}{\mu}\right)^{1/2}E\cdot \chi, \psi\cdot\chi\right\rangle
 = 0, \quad \psi \in H(\text{curl},\O), 
\nonumber \\
&&\qquad (\text{curl} ~ E, \varphi) + (\mu \frac{\p H_2}{\p t}, \varphi) = 0, \quad \varphi
\in L^2(\O).\label{eq_14}
\end{eqnarray}
Next, test \eqref{mod.1c} against $v^s\in [H^1(\O_p)]^2$ and \eqref{mod.1da} 
against  $v^f\in H(\text{div},\O_p)$ 
and use the integration by parts formula \eqref{by_parts2} 
and  the boundary condition \eqref{modf.6}. Setting 
\begin{eqnarray}\label{def_P}
\hskip4cm {\mathcal P} =\begin{pmatrix}
\rho_b I_d &\rho_f I_d \\
\rho_f I_d& m I_d  
\end{pmatrix}
\end{eqnarray}
where $I_d$ is the identity matrix in $R^d$ ($d=2$ here), we get the equation
\begin{eqnarray}
&& ({\mathcal P} \frac{\p^2 u}{\p t^2}, v)_{\O_p} 
+ ( \frac{\eta}{\kappa_0} \frac{\p u^f}{\p t}, v^f)_{\O_p} 
+ {\mathcal A}(u,v) + \left\langle {\mathcal D} S_{\G_p}(\frac{\p u}{\p t}),
S_{\G_p}(v)\right\rangle_{\G_p}\label{eq_15}\\
&&\hskip1cm   = (F,v)_{\O_p}, \quad v=(v^s, v^f)\in 
   [H^1(\O_p)]^2 \times H(\text{div},\O_p)\nonumber. 
\end{eqnarray}  

In \eqref{eq_15} $F= (F^s, F^f)$ and $A(u,v)$ is the bilinear form defined as 
\begin{eqnarray}
&&{\mathcal A}(u,v) = \sum_{l,m}\left(\tau_{lm}(u), \varepsilon_{lm}(v^s)\right)_{\O_p} - \left(p_f(u), 
\nabla\cdot v^f\right)_{\O_p} =\left({\bf M} ~ \tilde \epsilon(u),\tilde
\epsilon(v)\right)_{\O_p},\nonumber\\
&&\hskip2cm  \ u, v\in [H^1(\O_p)]^2\times
H(\text{div},\O_p),\label{def_A}
\end{eqnarray}
where the matrix ${\bf M}$  in \eqref{def_A} is given by  
\begin{eqnarray}\label{def_M}
\hskip4cm{\bf M} =\begin{pmatrix}
\lambda_c + 2 G&\lambda_c&\alpha~ K_{av}&0\\
\lambda_c&\lambda_c + 2 G& \alpha ~ K_{av} & 0\\
\alpha ~ K_{av}& \alpha ~ K_{av}& K_{av}&0\\
0&0&0&4 G
\end{pmatrix},
\end{eqnarray}
and  $\tilde\varepsilon(u) = \left(
\e_{11}(u^{(s)}), \e_{22}(u^{(s)}), \nabla\cdot u^{(f)}, 
\e_{12}(u^{(s)}\right)^t$.
 Furthermore, we assume that  ${\mathcal P}$ and ${\bf M}$ 
 are positive definite since they are associated with the kinetic and 
 strain energy densities, respectively. We also assume that the entries in 
 these two matrices are bounded below and above by positive constants.
 
Our continuous weak formulation is stated as follows: find 
$(E, H_2, u^s, u^f)\in H(\text{curl},\O)\times L^2(\O) 
\times [H^1(\O_p)]^2 \times H(\text{div},\O_p)$ satisfying  \eqref{eq_13},
\eqref{eq_14} and \eqref{eq_15}. 

Let us analyze the uniqueness of the solution of our continous problem. Set 
$F = 0$ in \eqref{eq_15} and  set to zero the initial conditions in
\eqref{icond}. Then 
choose $\varphi = H_2$ in \eqref{eq_14} to get 
\begin{eqnarray}
&&(\text{curl} ~ E, H_2) + (\mu \frac{\p H_2}{\p t}, H_2) = 0.\label{eq_16}
\end{eqnarray}
Also, choose $\psi = E$ in \eqref{eq_13} and use \eqref{eq_16} to obtain
\begin{eqnarray}
&& (\varepsilon \frac{ \p E}{\p t}, E)
 + (\sigma E, E) +  (\mu \frac{\p H_2}{\p t}, H_2) 
 + \left( L_0 \frac{\eta}{\kappa_0} u^f, E\right)_{\O_p}
 \label{eq_17}\\
&&\hskip3cm + \left\langle \left( \frac{\varepsilon}{\mu}\right)^{1/2}E\cdot \chi, E\cdot\chi\right\rangle
 = 0.\nonumber
\end{eqnarray} 
Next, choose $v^s = \frac{\p u^s}{\p t}, v^f = \frac{\p u^f}{\p t}$ in \eqref{eq_15}
and add the resulting  equation to  \eqref{eq_17} to get 

\begin{eqnarray}
&& \frac12 \frac{d}{dt}\left[ \left({\mathcal P} \frac{\p u}{\p t}, \frac{\p u}{\p
t}\right)_{\O_p} + \left( {\bf M} ~ \tilde \epsilon(u),\tilde
\epsilon(u)\right)_{\O_p} + (\varepsilon  E, E)  
+  (\mu  H_2, H_2)\right]\label{eq_22}\\
&&\qquad + \Phi(E,\frac{\p u^f}{\p t})  + (\sigma E, E)_{\O_a}
+ \left\langle \left( \frac{\varepsilon}{\mu}\right)^{1/2}E\cdot \chi, E\cdot\chi\right\rangle 
+  \left\langle {\mathcal D} S_{\G_p}(\frac{\p u}{\p t}),
S_{\G_p}(\frac{\p u}{\p t})\right\rangle_{\G_p} = 0,\nonumber 
\end{eqnarray} 
where 
\begin{eqnarray}
\Phi(E,\frac{\p u^f}{\p t}) = (\sigma E, E)_{\O_p} 
+ \left( L_0 \frac{\eta}{\kappa_0} u^f, E\right)_{\O_p} 
+ ( \frac{\eta}{\kappa_0} \frac{\p u^f}{\p t}, \frac{\p u^f}{\p
t})_{\O_p}.\label{def_Phi} 
\end{eqnarray}
Set  
\begin{eqnarray}
A_{\text{min},p} = \text{inf}\,\, \{A(x_1, x_3), (x_1,x_3)\in \O_p\}, \,\, 
A=\sigma, \kappa_0,\,\,   \kappa_{0,\text{max}} = \text{sup}\hskip.1cm  \{\kappa_0(x_1,
x_3), (x_1, x_3)\in \O_p\}.\nonumber    
\end{eqnarray}
Assume that
\begin{eqnarray}
&&\sigma_{\text{min},p}>0, \quad  \kappa_{0,\text{min},p}>0, \quad
\kappa_{0,\text{max},p}<\infty,\label{eq_20a}\\
&&C_1 = \text{min}\left(\sigma_{\text{min},p}  - 
\frac{L_0 \eta}{2 \kappa_{0,\text{max}},p}, 
\frac{\eta}{\kappa_{0,\text{max}},p}- 
\frac{L_0 \eta}{2 \kappa_{0,\text{max}},p}\right) > 0.\label{eq_20}
\end{eqnarray}
Thus, 
\begin{eqnarray}
\Phi(E,\frac{\p u^f}{\p t}) \ge C_1 \left(\|E\|^2_{0,\O_p} +
\|u^f\|^2_{0,\O_p}\right).\label{eq_21}
\end{eqnarray}
Next, since the matrix $\bf M$ is positive definite, Korn' second inequality 
\cite{duvaut-lions,nitsche81}  implies that 
\begin{eqnarray}
&&\left( {\bf M} ~ \tilde \epsilon(u),\tilde \epsilon(u)\right)_{\O_p} \ge C_2 
\left(\|u^s\|^2_{1,\O_p}\ + \|u^f\|^2_{H(\text{div},\O_p}\right) \label{eq_23}\\
&&\hskip3.2cm - C_3 \left(\|u^s\|^2_{0,\O_p}\ + \|u^f\|^2_{0,\O_p}\right).\nonumber
\end{eqnarray}
Set 
\begin{eqnarray*}
{\mathcal Z} = [H^1(\O_p)]^2 \times H(\text{div},\O_p),
\end{eqnarray*}
provided with the natural norm
\begin{eqnarray*}
\|u\|_{\mathcal Z} = \left( \|u^s\|^2_{1,\O_p} + \|u^f\|^2_{H(\text{div},\O_p}
\right).
\end{eqnarray*}
Then choose a constant $\zeta$ such that $\zeta > C_3$ and 
define the bilinear form
\begin{eqnarray*}
{\mathcal A}_{\zeta}(u,v) = {\mathcal A}(u,v) + \zeta (u,v),\label{def_a_zeta}
\end{eqnarray*}
so that  ${\mathcal A}_{\zeta}$ is ${\mathcal Z}$-coercive, {\it i.e.}, 
\begin{eqnarray}
{\mathcal A}_{\zeta}(u,u) \ge C_2 \|u\|^2_{\mathcal Z} ,\label{eq_25}
\end{eqnarray}
Next, add  to \eqref{eq_22} the inequality 
\begin{eqnarray*}
\frac{\zeta}{2} \frac{d}{dt} \|u\|^2_{0,\O_p} \le \frac{\zeta}{2}
\left(\|u\|^2_{0,\O_p} + \|\frac{\p u}{\p t}\|^2_{0,\O_p}\right),  
\end{eqnarray*}
integrate in time the resulting inequality  and use \eqref{eq_21} and
\eqref{eq_25} to obtain 
\begin{eqnarray}
&& \frac12\left[ \left({\mathcal P} \frac{\p u}{\p t}, \frac{\p u}{\p
t}\right)_{\O_p}(t) + C_2 \|u(t)\|^2_{\mathcal Z} + (\varepsilon  E, E)(t)  
+  (\mu  H_2, H_2)(t)\right]\label{eq_29}\\
&&\qquad + C_1 \int_0^t \left(\|E(s)\|^2_{0,\O_p} +
\|u^f(s)\|^2_{0,\O_p}\right) ds   
+ \int_0^t (\sigma E, E)_{\O_a}(s) ds \nonumber\\
&&\qquad + \int_0^t \left\langle \left(
\frac{\varepsilon}{\mu}\right)^{1/2}E\cdot \chi, 
 E\cdot\chi\right\rangle(s) ds   
+ \int_0^t \left\langle {\mathcal D} S_{\G_p}(\frac{\p u}{\p t}),
S_{\G_p}(\frac{\p u}{\p t})\right\rangle_{\G_p}(s) ds \nonumber\\
&&\qquad \le 
\int_0^t \left( \|u(s)\|^2_{0,\O_p} + \|\frac{\p u}{\p t}(s)\|^2_{0,\O_p} \right)
ds.\nonumber 
\end{eqnarray}   
Assume that $\varepsilon$ and $\mu$ are bounded below and above 
by positive constants, so
that
\begin{eqnarray}
0<\varepsilon_*\le \varepsilon\le \varepsilon^*<\infty, \quad 0<\mu_*\le \mu \le
\mu^*<\infty.\label{pos_epsilon_mu}
\end{eqnarray} 
Thus apply Gronwall's lemma in \eqref{eq_29},  note that all terms in the
left-hand side of \eqref{eq_29} are nonnegative and use \eqref{pos_epsilon_mu} 
to  conclude that
\begin{eqnarray*}
&&\|u(t)\|_{\mathcal Z} = 0, \quad  \|E(t)\|_{0}=0, \quad \|H_2(t)\|_{0}=0, \quad \forall t,
\end{eqnarray*} 
so that uniqueness holds for the solution of our differential problem.
Assuming sufficient regularity on the initial data and the external sources, 
 existence can be derived using the compacteness argument of Lions \cite{lions_69}
with an argument similar to that given in \cite{santos_86}. For brevity the
argument is ommited. 
The result is summarized in the following theorem.
\begin{theorem}\label{Theorem1}
%\label{thrm3.1}
Assume the validity of 
 \eqref{eq_20} and \eqref{pos_epsilon_mu} and that the matrices ${\bf M}$ and  
 ${\mathcal D}$ are positive definite. Then 
  there exists a unique  solution of problem \eqref{mod.1a}-\eqref{mod.1da} 
with the boundary conditions \eqref{modf.5}-\eqref{modf.6} and the initial
conditions \eqref{icond}.
\end{theorem}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%  SECTION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



 
\section{A continuous-time finite element method for the PSVTM-mode. Rectangular case.}
Let $\Tau^h(\O)$ be a nonoverlapping quasiregular 
 partition of $\O=\O_p\cup\O_a$ into rectangles $\O_j$ 
of diameter bounded by $h$ such that $\overline\Omega=\cup^J_{j=1}\overline \O_j$. 
Denote by $\xi_j$ and $\xi_{jk}$ the midpoints of $\G_j=\p \O_j\cap \G$ 
and $\G_{jk}= \G_{kj} = \p \O_j\cap\p \O_k $,
respectively. 
%Let us denote by $\nu_{jk}$ the unit outer normal on $\G_{jk}$ 
%from $\O_j$ to $\O_k$ and by $\nu_j$ the unit outer normal to $\G_j$.
%Let $\chi_j$ and $\chi_{jk}$ be unit tangents on $\G_j$ and 
%$\G_{jk}$ so that $\{\nu_j,\chi_j\}$ and 
%$\{\nu_{jk},\chi_{jk}\}$
%are orthonormal systems on $\G_j$ and $\G_{jk}$, respectively.


To approximate the electromagnetic fields $E, H_2$ we will employ 
the mixed finite element space    ${\mathcal V}^{h}\times {\mathcal W}^{h}$,  defined 
  as follows \cite{nedelec,pmonk_94}:
\begin{eqnarray*}
&&{\mathcal V}^{h}=\{\psi \in H(\text{curl},\O):\psi|\O_{j}\in {\mathcal V}^h_{j}
\equiv P_{0,1}(\O_j)\times P_{1,0}(\O_j)\},\\
&&{\mathcal W}^{h}=\{\varphi \in L^2(\O): \varphi|\O_{j}\in 
{\mathcal W}^h_{j} \equiv  P_0(\O_j)\}.
\end{eqnarray*}
Here, $P_{s,t}(\O_j)$   denote the polynomials of degree 
not greater than $s$ in $x$ and not greater than $t$ in $z$ on $\O_j$, while 
$P_0$ denote the constants on $\O_j$.
The functions in ${\mathcal V}^h$ have continuous tangential 
components across the internal boundaries $\G_{jk}$. Also, $\text{curl} {\mathcal V}^h\subset 
{\mathcal W}^{(h)}$ . 


Following \cite{pmonk_94}, the degrees of freedom for ${\mathcal V}^h$ 
are defined in the following way. 
Let $\O_{j}$
be a general element of the partition $\Tau^h(\O)$ and let
$\psi\in[H^1(\O_{j})]^2$. Then define the following moments on $\G_{jk}$:
\begin{eqnarray}\label{moments}
&&M_{\G_{jk}}(\psi)=\left\{\la\psi\cdot\tau,f\ra_{\G_{jk}};\quad f\in
P_{0}(\G_{jk})\right\}.
\end{eqnarray}
Note that  \eqref{moments}  are
curl--conforming and unisolvent for elements in ${\mathcal V}^h$.

To approximate each component of  the solid displacement vector
we employ the nonconforming finite element space  ${\mathcal NC}^h$  
as in \cite{dssy-ncell}, while to approximate the fluid displacement  
vector we choose ${\mathcal M}^{h}$, the vector part of the 
Raviart-Thomas-Nedelec space 
\cite{rath,nedelec} of zero order. 
More precisely, set 
$$
\widehat R=[-1,1]^2, \qquad
\widehat {\mathcal NC}(\widehat R)=\mbox{Span}\{1,\hat x, \hat z,\widetilde  \alpha(\hat x)-
\widetilde \alpha(\hat z)\},\quad
\widetilde \alpha(\hat x) =  \hat x^2 - \frac{5}{3}  \hat x^4.
$$
with the degrees of freedom being the values at the 
midpoint of each edge of $\widehat R$. 
Next, for each $\O_j\subset\O_p$, let
$F_{\O_j} : \widehat R\rightarrow \O_j$
be an invertible affine mapping such that $F_{\O_j}(\widehat R)=\O_j$, and define
\begin{eqnarray*}
{\mathcal NC}^h_j &=&\{v:\, v = \widehat v\circ F_{\O_j}^{-1},\, 
\widehat v\in \widehat {\mathcal NC}(\widehat R)\}.\nonumber
\end{eqnarray*}
Thus,  
\begin{eqnarray*}
&&{\mathcal NC}^{h} = \{ v \, :  v_j=  v|_{\O_j}\,  \in 
{\mathcal NC}^h_j, ~ 
v_j(\xi_{jk}) =v_k(\xi_{jk})\,\,  \forall (j,k)\} ,\nonumber\\
&&{\mathcal M}^{h}=\{w \in H(\text{div},\O_p):w|\O_{j}\in {\mathcal M}^h_{j}
\equiv P_{1,0}(\O_j)\times P_{0,1}(\O_j)\}.\nonumber
\end{eqnarray*}
To state the approximating properties of the finite 
element spaces defined above we introduce the following 
four projection operators.

First, let 
\[
[H^1_h(\O)]^2=\{\psi:\psi|\O_{j}\in[H^1(\O_{j})]^2\},
\]
with $[H^1_h(\O_p)]^2$ defined in similar fashion  and  
 if $\G_{jk,p}$ denotes any inner interface $\G_{jk}$ in $\O_p$  let 
\begin{eqnarray*} 
&&\tilde \Lam^h = \left\{\tilde \lam^h : \tilde \lam^h_{{jk}} = \text{tr}_{~\G_{jk,p}}(\tilde
\lam^h|_{\O_j})\in
[P_0(\G_{jk,p} )]^2\equiv\tilde \Lam^h_{jk},
\quad\tilde\lam^h_{jk} + \tilde\lam^h_{kj} = 0 \right\},\nonumber
\end{eqnarray*}
where  $P_0(\G_{jk,p})$ denotes the constant  functions defined on
$\G_{jk,p}$.

\noindent
{\emph Remark}: Note that there are two     
copies of $[P_0(\G_{jk,p} )]^2$ assigned 
to each $\G_{jk,p}$, one from $\O_j$ to $\O_k$ and another  
from $\O_k$ to $\O_j$.

Then  we define the  projections
\begin{eqnarray}
&&\hskip-.9cm \Pi_h:H(\text{curl},\O)\cap[H^1_h(\O)]^2\to {\mathcal V}^h:   
\la(\psi-\Pi_h\psi)\cdot\chi, 1\ra_{B}=0,\, B=\G_{jk} \, \mbox{ or } \, \G_j,\label{def_pih}\\
&&\hskip-.9cm P_h : L^2(\O)\rightarrow  {\mathcal W}^{h}: 
\hskip3cm (P_h w - w, \varphi)=0, \,    
\varphi \in {\mathcal W}^{h}, \label{def_ph}\\ 
&&\hskip-.9cm R_h: [H^2(\O_p)]^2\rightarrow [\NC^h]^2:
\hskip2cm (v_i^{s} - R_h v_i^{s})(\xi) = 0, \, \xi = \xi_{jk}\mbox{ or } \xi_j,\label{def_rh}\\
&&\hskip6.5cm \text{for}\quad   v^s=(v^s_1, v^s_2),\nonumber\\ 
&&\hskip-.9cm Q_h: [H^1(\O_p)]^2\rightarrow {\mathcal M}^h:
\hskip1cm \la (v^{f}-Q_h v^{f})\cdot\nu, 1\ra_{B}=0;  
\, B=\G_{jk,p} \, \mbox{ or } \, \G_j,\label{def_qh}\\
&&\hskip-.9cm S_h : [H^2(\O_p)]^2\times H^1(\text{div};\O_p) \rightarrow \tilde \Lam^h:\qquad
\left\langle \tau(v)\nu - S_h(v), 1\right\rangle_{B} = 0,\label{def_sh}\\
&&\hskip8cm v=(v^{s}, v^{f}),  B=\G_{jk,p} \, \mbox{ or } \, \G_j.\nonumber
\end{eqnarray}
Let us define the  broken norms
\[
\|v\|^2_{s,h,\O_p} = \sum_{\O_{j}\subset \O_p} \|v\|^2_{s,\O_j}. 
\]
The   approximation properties of these operators can be 
stated as follows \cite{nedelec,santos_sinum_07}:
\begin{eqnarray}
&&\|\psi-\Pi_h\psi\|_0\leq C h\|\psi\|_{1},\,\psi\in [H^1(\O)]^2,\label{aprox_pih1}\\
&&\|\text{curl}(\psi-\Pi_h\psi)\|_0\leq C h\|\text{curl}~\psi\|_{1},\,
\psi\in [H^1(\O)]^2,\, \text{curl}\psi \in H^1(\O),\label{aprox_pih2} \\
&&\|P_h \var- \var\|_0\leq C h\|\var\|_1, \qquad \forall \var\in H^1(\O),\label{aprox_ph}\\
&&\left[\|v^s- R_h v^{s}\|_{\O_p}
+ h \|v^{s}-R_h v^{s}\|^2_{1,h,\O_p}
+ h^2 \|v^{s}
- R_h v^{s}\|^2_{2,h,\O_p} 
\nonumber\right.\\
&&\qquad
+  h^{\frac12} \bigg(\sum_{\O_{j}\subset \O_p}\|v^{s}-R_h v^{s}\|^2_{0,\p \O_j}
\bigg)^{\frac12}
+ h^{\frac{3}{2}}\bigg(\sum_{\O_{j}\subset \O_p}\|\tau(v_j)\nu_j -
S_h v_j\|^2_{0,\p \O_j}\bigg)^{1/2}\label{approx_rh} \\
&&\quad
\le C  h^2 \left(\|v^{s}\|_{2,\O_p}  + 
\|\nabla\cdot v^{f}\|_{1,\O_p}\right),\,v=\left(v^{s}, v^{f}\right)
\in [H^2(\O_p)]^2\times H^1(\text{div}, \O_p)\nonumber,\\
&&\|Q_h v^f - v^f\|_{0,\O_p}\leq C h\|\v^f\|_{1,\O_p},\,v^f\in [H^1(\O_p)]^2,\label{aprox_qh1}\\
&&\|\nabla\cdot(v^f- Q_h v^f)\|_{0,\O_p}\leq C h\|\nabla\cdot v^f\|_{1,\O_p},\,v^f 
\in H^1(\text{div},\O_p)\label{aprox_qh2}.
\end{eqnarray}
Note that since $\text{curl} ~\psi\in {\mathcal W}^h\, \forall \psi\in {\mathcal V}^h$,
  it follows from \eqref{def_ph} that 
\begin{eqnarray}\label{orthog_1}
&&\left(P_h f - f, \text{curl} ~\psi \right) = 0, \quad \forall \psi \in {\mathcal V}^h.
\end{eqnarray}
Also note the orthogonality  property for functions on ${\mathcal NC}^h$:
\begin{eqnarray}\label{orth}
\< v^s_j -v^s_k, 1\>_{\G_{jk}} = 0 \mbox{ for  
all interior  interfaces } \G_{jk} ,
\quad v^s \in \NC^h.
\end{eqnarray}


Set 
\begin{eqnarray}\label{bilinear3}
&&{\mathcal A}_h(u,v)=\sum_{\O_{j}\subset \O_p}\left[\sum_{l,m}\left(\tau_{lm}(u),
\e_{lm}(v^{(s)})\right)_{\O_j} 
- \left(p_f(u),\nabla\cdot v^f)\right)_{\O_j}\right]\\
&&\qquad\qquad =\sum_{\O_{j}\subset \O_p} \left({\bf M} ~ \tilde \epsilon(u), 
\tilde \epsilon(v)\right)_{\O_j},  \nonumber
\end{eqnarray}
and
\begin{eqnarray}
&&\Theta_h\left( (E,H_2,u^s,u^f), (\psi,\varphi,v^s,v^f)\right) 
=(\varepsilon  \frac{\p E}{\p t}, \psi) +
 (\sigma E^, \psi) -(H_2, \text{curl} \psi) \label{def_Thetah}\\
&& \qquad + \left( L_0 \frac{\eta}{\kappa_0} \frac{\p  u^f}{\p t}, \psi\right)_{\O_p} 
 +(\text{curl} E, \varphi) + \left\langle \left(\frac{\varepsilon}{\mu}\right)^{1/2}  E\cdot \chi, 
\psi\cdot\chi\right\rangle \nonumber\\
&&\qquad  + (\mu \frac{\p H_2}{\p t}, \varphi)
+ \left({\mathcal P} \frac{\p^2 u}{\p t^2} , v\right)_{\O_p} + \left( \frac{\eta}{\kappa_0} \frac{\p 
u^f}{\p t}, v^f\right)_{\O_p}
+  {\mathcal A}_h(u,v) + \left\langle {\mathcal D} S_{\G}(\frac{\p u}{\p t}),
 S_{\G}(v)\right\rangle_{\G_p}.\nonumber
\end{eqnarray}



Let $J= (0,T)$ and 
\[
{\mathcal Y}^h = {\mathcal V}^h \times {\mathcal W}^h \times 
({\mathcal NC}^h)^2 \times {\mathcal M}^h.
\]   
Then the {\it global} Galerkin procedure is defined as follows: 
find $\left(E^{h}, H^{h}_2, u^{s,h}, u^{f,h}\right) :J \rightarrow  
{\mathcal Y}^h$ 
such that 
\begin{eqnarray}    
&&\Theta_h\left( (E^{h},H^{h}_2,u^{s,h},u^{f,h}), (\psi,\varphi,v^s,v^f)\right)
 = (F^{s}, v^s)_{\O_p} + (F^{f}, v^f)_{\O_p}\label{eq_33a}\\
&&\hskip4cm  (\psi, \varphi, v^s, v^f)\in {\mathcal Y}^h,\nonumber\\
&&\text{with the initial conditions}\nonumber\\
&&E^{h}(t=0)\approx E_0,\quad  H^{h}_2(t=0) \approx H_{2,0}\label{eq_33b}\\ 
&&u^{s,h}(t=0)\approx u^s_0,\quad u^{f,h}(t=0)\approx u^f_0.\label{eq_33c}
\end{eqnarray}
To analyze the uniqueness for \eqref{eq_33a}-\eqref{eq_33c}, let us introduce the space 
\[
{\mathcal Z}^h = [ H^1_h(\O_p)]^2 \times H(\text{div}, \O_p) 
\]
provided with the norm
\[
\|u\|^2_{{\mathcal Z}^h} = \left( \|u^s\|^2_{1,h,\O_p} 
+ \|u^f\|^2_{H(\text{div}, \O_p)}\right)^{1/2}.
\]
As in the continuous case, let us define  the bilinear form  
\begin{eqnarray}
{\mathcal A}_{\zeta,h}(u,v) = {\mathcal A}_{h}(u,v) + \zeta (u, v),\quad u,v \in {\mathcal
Z}^h, 
\end{eqnarray}
with $\zeta > C_3$, so that
\begin{eqnarray}
&&{\mathcal A}_{\zeta,h}(u,u)\ge C_2 \|u\|^2_{{\mathcal
Z}^h},\label{coercivity_ah}\\
&&{\mathcal A}_{\zeta,h}(u,v)\le C \|u\|_{{\mathcal Z}^h} \|v\|_{{\mathcal
Z}^h}.\label{continuity_ah}
\end{eqnarray} 
Then if $u^h = (u^{s,h}, u^{f,h})$, 
a repetition of the argument leading
to \eqref{eq_29} implies that 
\begin{eqnarray*}
&&\|u^h\|_{{\mathcal Z}^h} = 0, \quad  \|E^h(t)\|_{0}=0, \quad \|H^h_2(t)\|_{0}=0, 
\quad \forall t,
\end{eqnarray*} 
so that uniqueness holds for \eqref{eq_33a}-\eqref{eq_33c}. Existence follows from
finite-dimensionality. Thus we can state the following theorem.
\begin{theorem}\label{Theorem2}
%\label{thrm3.1}
Under the hypothesis  of the positive definitess of the matrices 
${\bf M}$and  ${\mathcal D}$  and the validity of \eqref{eq_20} and
\eqref{pos_epsilon_mu},  there exists a unique  solution  of \eqref{eq_33a} for every
choice of the initial conditions \eqref{eq_33b}and  \eqref{eq_33c}.
\end{theorem}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%  SECTION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



 
\section{Apriori error estimates for the PSVTM-mode}

Now we derive the  error estimate for the procedure \eqref{eq_33a} stated in the
following theorem.
\begin{theorem}\label{Theorem3}
%\label{thrm3.1}

Assume that  the matrices 
${\bf M}$, ${\mathcal D}$ are positive definite and  the validity of \eqref{eq_20} and
\eqref{pos_epsilon_mu}. Also assume that $\varepsilon$ and $\mu$ are piecewise constant
on $\G$. Then   
the following {\it a priori}  error estimate holds: 
\begin{eqnarray}
&& \|E -E^{h}\|_{L^{\infty}(J,L^2(\O))} +   \|H_2 -H_2^{h}\|_{L^{\infty}(J,L^2(\O))} 
+ \|u^s -u^{s,h}\|_{L^{\infty}(J,H^1_h(\O_p))}\nonumber\\ 
&&\quad + \|u^f -u^{f,h}\|_{L^{\infty}(J,H(\text{div},\O_p))} 
+\| \frac{\p (u^s-u^{s,h})}{\p t})\|_{L^{\infty}(J,L^2(\O_p))} \label{apriori}\\
&&\quad+\| \frac{\p (u^f-u^{f,h})}{\p t})\|_{L^{\infty}(J,L^2(\O_p))} 
   +  \|u^s -u^{s,h}\|_{L^{2}(J,L^2(\G_p))} \nonumber\\
&&\qquad +  \|(u^f-u^{f,h})\cdot \nu \|_{L^{2}(J,L^2(\G_p))} \nonumber\\
&&\qquad\qquad\le C \left[ h\left(N_0 +  N_1\right)
+  h^{1/2}\left( M_0 +  M_1 \right)\right],\nonumber
\end{eqnarray}
where
\begin{eqnarray*}
&&M_1^2 = \|u^s\|^2_{L^{\infty}(J,H^2(\O_p))} 
+ \|\frac{\p u^s}{\p t}\|^2_{L^{\infty}(J,H^2(\O_p))} 
+ \|u^f\|^2_{L^{\infty}(J,H^1(\text{div},\O_p))}\\ 
&&\qquad + \|\frac{\p u^f}{\p t}\|^2_{L^2(J,H^1(\text{div},\O_p))},\\
&&N_1^2 = \|E\|^2_{L^{\infty}(J,H^1(\O))} 
+ \|\text{curl} ~ E\|^2_{L^{\infty}(J,H^1(\O))} 
+ \|\frac{\p E}{\p t}\|^2_{L^{\infty}(J,H^1(\O))} 
+\|\frac{\p H_2}{\p t}\|^2_{L^{2}(J,H^1(\O))},\nonumber\\
&& M_0 = \|u^s(0)\|_{2,\O_p} + \|u^f(0)\|_{H^1(\text{div},\O_p)} +
 \|\frac{\p u^s}{\p t}(0)\|_{2,\O_p} 
 + \|\frac{\p u^f}{\p t}(0)\|_{H^1(\text{div},\O_p},\\
&& N_0 = \|E(0)\|_{1} +\|H_2(0)\|_{1}. 
\end{eqnarray*}
 

\end{theorem}
\begin{proof}
Set  
\begin{eqnarray*}\label{error1}
&& \delta = \left(E^s -E^{h}, H_2 - H^{h}_2, u^s - u^{s,h}, 
u^f - u^{f,h}\right)\equiv\left(\delta^E,\delta^H_2, \delta^s, \delta^f \right),
\\
&&\gamma = \left (\Pi_h E- E^{h}, P_h H_2 -H^{h}_2 ,
 {\mathcal  R}_h u^s - u^{s,h}, Q_h u^f - u^{f,h}\right) 
 \equiv\left(\gamma^E, \gamma^H_2, \gamma^s, \gamma^f \right),\\ 
&&\text{and}\\
&&\delta^B = (\delta^s, \delta^f), \quad  \gamma^B = (\gamma^s, \gamma^f).
\end{eqnarray*}
First, use integration by parts element by element and the boundary conditions
\eqref{modf.5}-\eqref{modf.6}  to see that, for $\psi\in H(\text{curl},\O)$,
$\varphi\in L^2(\O)$ and  $v=(v^s,v^f)\in \left[L^2(\O_p)\right]^4$ such that 
$v^s|_{\O_j}\in [H^1(\O_j)]^2$ and  $ v^f|_{\O_j}\in H(\text{div}, \O_j)$,  
\begin{eqnarray}
&&\Theta_h\left( (E, H_2, u^s, u^f), (\psi,\varphi,v^s,v^f)\right) = 
  (F^{s}, v^s)_{\O_p} + (F^{f}, v^f)_{\O_p} \label{eq_42}\\
&&\qquad + \sum_{\O_{j}\subset \O_p}\left[\la\tau(u)\nu, v^s\ra_{\p\O_j\setminus \G_p} -
 \la p_f(u), v^f\cdot\nu\ra_{\p\O_j\setminus\G_p}\right].\nonumber
\end{eqnarray}
Now subtract \eqref{eq_42} from \eqref{eq_33a} to 
see that, for $\left(\psi,\varphi,v^s,v^f\right)\in {\mathcal Y}^h$  
\begin{eqnarray}
&&\Theta_h\left( (\delta^E,\delta^H_y, \delta^s, \delta^f), (\psi,\varphi,v^s,v^f)\right) =  
\sum_{\O_{j}\subset \O_p}\left[\la\tau(u)\nu, 
 v^s \ra_{\p\O_j \setminus \G_p}\label{eq_44} \right.\\
&&\left. \hskip7cm - 
 \la p_f(u), v^f\cdot\nu\ra_{\p\O_j\setminus\G_p}
\right].\nonumber
\end{eqnarray}
Next, by  hypothesis $u^s\in [H^2(\O_p)]^2, ~ u^f\in H^1(\text{div},\O_p)$ and 
consequently  $p_f(u)\in H^{1/2}(\p\O_j)$. Also, $v^f\in {\mathcal M}^h\subset
H(\text{div},\O)$ and consequently  $v^f\cdot\nu$ 
is continuous across  the interior interfaces of the elements $\O_j$. Thus, 
\begin{eqnarray}\label{eq_43}
\sum_{\O_{j}\subset \O_p} \la p_f(u), v^f\cdot \nu \ra_{\p\O_j \sm\G_p} = 0.
\end{eqnarray}
Also  use that $v^s\in \NC^h$ is  constant on $\G_{jk}$ and that  $S_h(u)$ changes sign on each side of $\G_{jk}$ (c.f. \eqref{def_sh}) to see that 
\begin{eqnarray}\label{eq_45}
\sum_{\O_{j}\subset \O_p} \la S_h(u), v^s \ra_{\p\O_j\sm\G_p} = 0.
\end{eqnarray}
Thus, \eqref{eq_42} becomes
\begin{eqnarray}
&&\Theta_h\left((\delta^E,\delta^H_y, \delta^s, \delta^f) , (\psi,\varphi,v^s,v^f)\right) 
=\sum_{\O_{j}\subset \O_p} 
\la\tau(u)\nu - S_h(u), v^s\ra_{\p\O_j \setminus \G_p},\label{eq_46} \\
&&\hskip7cm \left(\psi,\varphi,v^s,v^f\right)\in {\mathcal Y}^h.\nonumber 
\end{eqnarray}
Now \eqref{eq_46} yields
\begin{eqnarray}
&&\Theta_h\left((\gamma^E,\gamma^H_2, \gamma^s, \gamma^f) , (\psi,\varphi,v^s,v^f)\right)\nonumber \\ 
&&\qquad = \Theta_h\left(\left(\Pi_h E - E, P_h H_2 - H_2, R_h u^s - u^s, 
Q_h u^f - u^f\right) , \left(\psi,\varphi,v^s,v^f\right)\right) \nonumber\\ 
&&\qquad\qquad  + \sum_{\O_{j}\subset \O_p} 
\la\tau(u)\nu - S_h(u), v^s\ra_{\p\O_j \setminus \G_p}, \quad
\left(\psi,\varphi,v^s,v^f\right)\in {\mathcal Y}^h\label{eq_47}
\end{eqnarray}
Next, choose $\psi =0, \varphi=0, v^s = \frac{\p \gamma^s}{\p t}, 
v^f = \frac{\p \gamma^f}{\p t}$ in \eqref{eq_47} and add to the resulting
equation the inequality
\begin{eqnarray*}
\frac{\zeta}{2} \frac{d}{dt} \|\gamma^B\|^2_{0,\O_p} \le \frac{\zeta}{2}
\left(\|\gamma^B\|^2_{0,\O_p} + \|\frac{\p \gamma^B}{\p t}\|^2_{0,\O_p}\right),  
\end{eqnarray*}
to obtain the inequality
\begin{eqnarray}
&&\frac12 \frac{d}{dt} \left[\left({\mathcal P} \frac{\p \gamma^B}{\p t}, 
\frac{\p \gamma^B}{\p t}\right)_{\O_p}  
+  {\mathcal A}_h(\gamma^B,\gamma^B)\right]
+ \left( \frac{\eta}{\kappa_0} \frac{\p \gamma^f}{\p t}, \frac{\p \gamma^f}{\p t}  \right)_{\O_p}
 + \left\langle {\mathcal D} S_{\G_p}(\frac{\p \gamma^B}{\p t}),
 S_{\G_p}(\frac{\p \gamma^B}{\p t})\right\rangle_{\G_p}\nonumber\\
&&\qquad \le \frac{\zeta}{2}
\left(\|\gamma^B\|^2_{0,\O_p} + \|\frac{\p \gamma^B}{\p t}\|^2_{0,\O_p}\right)
\nonumber\\
&&\qquad + \left({\mathcal P} (R_h u^s- u^s, Q_h u^f - u^f), 
\frac{\p \gamma^B}{\p t}\right)_{\O_p}  
+  {\mathcal A}_h((R_h u^s- u^s, Q_h u^f - u^f) , \frac{\p \gamma^B}{\p
t})\label{eq_48}\\
&&\qquad + \left( \frac{\eta}{\kappa_0} \frac{\p (Q_h u^f - u^f)}{\p t}, 
\frac{\p \gamma^f}{\p t}  \right)_{\O_p}
 + \left\langle {\mathcal D} S_{\G_p}((R_h u^s- u^s, Q_h u^f - u^f)),
 S_{\G_p}(\frac{\p \gamma^B}{\p t})\right\rangle_{\G_p}\nonumber\\
&& \qquad + \sum_{\O_{j}\subset \O_p} \la\tau(u)\nu - S_h(u), 
\frac{\p \gamma^s}{\p t}\ra_{\p\O_j \setminus \G_p},\nonumber 
\end{eqnarray}
which will be integrated in time 
 from  $0$ to $t$, but first we will bound the time
integral of the last five terms in the right-hand side of \eqref{eq_48}
First, using the approximation properties \eqref{approx_rh} and  \eqref{aprox_qh1}
\begin{eqnarray} 
&&\left|\left({\mathcal P} (R_h u^s- u^s, Q_h u^f - u^f),
\frac{\p \gamma^B}{\p t}\right)_{\O_p}\right| \le 
C \left(h^4 \|u^s\|^2_{2,\O_p} + h^2 \|u^f\|^2_{1,\O_p}\right.\nonumber\\ 
&&\hskip7.5cm \left. 
+ \|\frac{\p \gamma^B}{\p t}\|^2_{2,\O_p}\right),\nonumber 
\end{eqnarray}
so that 
\begin{eqnarray} 
&&\int_0^t \left|\left({\mathcal P} (R_h u^s- u^s, Q_h u^f - u^f),
\frac{\p \gamma^B}{\p t}\right)_{\O_p}(s)\right| ds\label{eq_63}\\
&&\qquad \le C \left[\int_0^t \|\frac{\p \gamma^B}{\p t}(s)\|^2_{0,\O_p} ds 
+ h^2 M_1 \right].\nonumber  
\end{eqnarray} 
Next,  using again \eqref{aprox_qh1}
\begin{eqnarray} 
&&\int_0^t \left|\left(\frac{\eta}{\kappa_)} 
\frac{\p Q_h u^f - u^f)}{\p t}, \frac{\p \gamma^f}{\p t}\right)_{\O_p}(s) ds\right|
\le \widehat \epsilon  \int_0^t \left(\frac{\p \gamma^f}{\p t},
\frac{\p \gamma^f}{\p t}\right)_{\O_p})(s) ds \label{eq_64}\\
&&\hskip7cm +  C  h^2 M_1.\nonumber  
\end{eqnarray} 
Also, using integration by parts in time, 
\begin{eqnarray}
&&\int_0^t {\mathcal A}_h((R_h u^s- u^s, Q_h u^f - u^f) , \frac{\p \gamma^B}{\p
t})(s) ds \label{eq_65}\\
&&\qquad = {\mathcal A}_h((R_h u^s- u^s, Q_h u^f - u^f) ,
\gamma^B)|_0^t\nonumber\\
&&\qquad - \int_0^t {\mathcal A}_h((\frac{\p R_h u^s- u^s}{\p t}, 
\frac{\p Q_h u^f - u^f}{\p t}) , \gamma^B)(s) ds.\nonumber
\end{eqnarray}
Let us bound each term in the right-hand side of \eqref{eq_65}. First, 
\begin{eqnarray}
&&\left|{\mathcal A}_h((R_h u^s- u^s, Q_h u^f - u^f) ,
\gamma^B)(0)\right| \label{eq_66}\\
&&\le C\left( \|(R_h u^s-
u^s)(0)\|_{1,h,\O_p}+\|(Q_h u^f- u^f)(0)\|_{H(\text{div},\O_p}\right) 
\|\gamma^B(0)\|_{{\mathcal Z}^h}\nonumber\\
&&\qquad \le C\left( \|\gamma^B(0)\|^2_{{\mathcal Z}^h} + h^2 M_0^2\right).\nonumber
\end{eqnarray}
Now choose the initial condition $u^{s,h}(0)\in {\mathcal NC}^h, u^{f,h}(0) \in
{\mathcal M}^h$ such that 
\begin{eqnarray}
A_{\zeta,h}(\delta^B(0),v) = 0, \quad v=(v^s,v^f)\in {\mathcal NC}^h \times
{\mathcal M}^h.\label{def_initial1}
\end{eqnarray}
Then choose $v = \delta^B(0) + (R_h u^s(0), Q_h u^f(0)) - (u^s(0), u^f(0))$
in \eqref{def_initial1} and use  the coercivity of  $A_{\zeta,h}$ in ${\mathcal
Z}^h$ to see that
\begin{eqnarray}
\|\delta^B(0)\|_{{\mathcal Z}^h }\le C h ~M_0, 
\end{eqnarray}
so that using the triangle inequality  and 
\eqref{approx_rh} and  \eqref{aprox_qh1} we obtain the bound
\begin{eqnarray}
\|\gamma^B(0)\|_{{\mathcal Z}^h }\le C h ~ M_0. 
\end{eqnarray}
Thus, \eqref{eq_66} becomes
\begin{eqnarray}
&&\left|{\mathcal A}_h((R_h u^s- u^s, Q_h u^f - u^f) ,
\gamma^B)(0)\right| \le C h^2 M_0^2.\label{eq_66a}
\end{eqnarray}
Next, 
\begin{eqnarray}
&&\left|{\mathcal A}_h((R_h u^s- u^s, Q_h u^f - u^f) ,
\gamma^B)|(t)\right| \label{eq_67}\\
&&\qquad \le \widehat \epsilon \|\gamma^B(t)\|^2_{{\mathcal Z}^h} +C  h^2 \left(
 \|u^s(t)\|^2_{2,\O_p} 
 + \|u^f(t)\|^2_{H^1(\text{div},\O_p)}\right),\nonumber
\end{eqnarray}
and 
\begin{eqnarray}
&&\int_0^t \left|{\mathcal A}_h((\frac{\p R_h u^s- u^s}{\p t}, 
\frac{\p Q_h u^f - u^f}{\p t}) , \gamma^B)(s) ds\right| \label{eq_68}\\
&&\qquad \le C\left(\int_0^t \|\gamma^B(s)\|_{{\mathcal Z}^h} ds  + h^2
M_1^2\right).\nonumber
\end{eqnarray}
Thus, collecting the estimates \eqref{eq_66a},\eqref{eq_67}and  \eqref{eq_68} we
get  the bound
\begin{eqnarray}
&&\int_0^t \left|{\mathcal A}_h((R_h u^s- u^s, Q_h u^f - u^f) , \frac{\p \gamma^B}{\p
t})(s) ds\right| \label{eq_69}\\
&&\qquad \le \widehat \epsilon \|\gamma^B(t)\|^2_{{\mathcal Z}^h} 
+ C\left(\int_0^t \|\gamma^B(s)\|_{{\mathcal Z}^h} ds  + h^2 (M_0^2 +
M_1^2)\right).\nonumber
\end{eqnarray}
Note that the next term in the right-hand side of \eqref{eq_48} can 
be bounded as follows:
\begin{eqnarray*}
&&\left|\left\langle {\mathcal D} S_{\G_p}((R_h u^s- u^s, Q_h u^f - u^f)),
 S_{\G_p}(\frac{\p \gamma^B}{\p t})\right\rangle_{\G_p}\right|\\
&&\qquad \le C\sum_j h^{1/2}\left(\|u^s\|_{1/2,\G_j} + \|u^f\cdot \nu_j\|_{1/2,\G_j} \right)
\left(\|\frac{\p \gamma^s}{\p t}\|_{0,\G_j} 
+ \|\frac{\p \gamma^f\cdot \nu_j}{\p t}\|_{0,\G_j}\right)\nonumber\\
&&\qquad \le \widehat \epsilon \left\langle {\mathcal D} S_{\G_p}(\frac{\p
\gamma^B}{\p t}), S_{\G_p}(\frac{\p \gamma^B}{\p t})\right\rangle_{\G_p} +
C h \left( \|u^s\|^2_{1,\O_p} + \|u^f\|^2_{1,\O_p}\right),\nonumber
\end{eqnarray*}
so that
\begin{eqnarray}  
&&\int_0^t \left|\left\langle {\mathcal D} S_{\G_p}((R_h u^s- u^s, Q_h u^f - u^f)),
 S_{\G_p}(\frac{\p \gamma^B}{\p t})\right\rangle_{\G_p}(s) ds
 \right|\label{eq_62}\\
&&\qquad \le \widehat \epsilon \int_0^t \left\langle {\mathcal D} S_{\G_p}(\frac{\p
\gamma^B}{\p t}), S_{\G_p}(\frac{\p \gamma^B}{\p t})\right\rangle_{\G_p}(s) ds +
C h ~M_1^2.\nonumber 
\end{eqnarray}
Next, using integration by parts in time in the last term in the right-hand side
of \eqref{eq_48},
\begin{eqnarray}
&&\int_0^t \la\tau(u)\nu - S_h(u), 
\frac{\p \gamma^s}{\p t}\ra_{\p\O_j \setminus \G_p}(s) ds =
\la\tau(u)\nu - S_h(u), 
\gamma^s\ra_{\p\O_j \setminus \G_p}|_0^t \nonumber\\
&&\hskip3cm  - 
\int_0^t \la\tau(\frac{\p u}{\p t})\nu - S_h(\frac{\p u}{\p t}), 
\gamma^s\ra_{\p\O_j \setminus \G_p}(s) ds.\label{eq_49}
\end{eqnarray}
Next, since $(\tau(u)\nu - S_h(u))(s)$ has average value zero 
on $\p\O_j \setminus \G_p$, if  $q^s(s)$ is the average value of $\gamma^s(s)$ 
on $\O_j$ for $s=0,\cdots,t$, using 
the trace inequality and Korn's  second inequality, \cite{duvaut-lions,nitsche81},
\begin{eqnarray}
&& \left|\la\tau(u)\nu - S_h(u), \gamma^s\ra_{\p\O_j \setminus \G_p}\right(0)|\nonumber\\
&&\qquad =\left|\la\tau(u)\nu - S_h(u), \gamma^s- q^s\ra_{\p\O_j \setminus
\G_p}\right(0)|\nonumber\\
&&\qquad \le  \|(\tau(u)\nu - S_h(u))(0)\|_{0,\p\O_j} 
\|(\gamma^s- q^s)(0)\|_{0,\p\O_j}\label{eq_55}\\
&&\qquad \le C h^{1/2} \left(\|u^s(0)\|_{2,\O_j}  + \|\nabla\cdot u^f(0)\|_{1,\O_j} \right)
\|(\gamma^s- q^s)(0)\|^{1/2}_{0,\O_j} \|(\gamma^s- q^s)(0)\|^{1/2}_{1,\O_j}\nonumber\\
&&\qquad \le  C h  \left(\|u^s(0)\|_{2,\O_j}  + \|\nabla\cdot u^f(0)\|_{1,\O_j} \right) 
\|\gamma^s(0)\|_{1,\O_j}.
\end{eqnarray}
Thus, adding over $j$ in \eqref{eq_55} we obtain
\begin{eqnarray}
&& \sum_j \left|\la\tau(u)\nu - S_h(u), \gamma^s\ra_{\p\O_j \setminus
\G_p}(0)\right|
\le  C h^2  ( M_0^2 + M_1^2).\label{eq_56}
\end{eqnarray}
Similarly, 
\begin{eqnarray}
&&\left|\la\tau(u)\nu - S_h(u), \gamma^s-
q^s\ra_{\p\O_j\setminus\G_p}(t)\right|\label{eq_57}\\
&&\qquad \le  \widehat\epsilon \|\gamma^s(t)\|^2_{1,\O_j} + 
C h^2  \left(\|u^s(t)\|^2_{2,\O_j}  + \|\nabla\cdot
u^f(t)\|^2_{1,\O_j}\right),\nonumber
\end{eqnarray}
so that adding over $j$ in \eqref{eq_57} we get
\begin{eqnarray}
&&\sum_j\left|\la\tau(u)\nu - S_h(u), \gamma^s-
q^s\ra_{\p\O_j\setminus\G_p}(t)\right|\label{eq_58}\\
&&\qquad \le  \widehat\epsilon \|\gamma^s(t)\|^2_{1,h} + 
C h^2 M_1^2.\nonumber
\end{eqnarray}
In a similar fashion,
\begin{eqnarray}
&&\sum_j \left|\int_0^t \la\tau(\frac{\p u}{\p t})\nu - S_h(\frac{\p u}{\p t}), 
\gamma^s\ra_{\p\O_j \setminus \G_p}(s) ds\right|\label{eq_59}\\
&&\qquad \le C\left( \int_0^t \|\gamma^s(s)\|^2_{1,h} ds + h^2
M_1^2\right).\nonumber
\end{eqnarray}
Combining the estimates \eqref{eq_56}, \eqref{eq_58} and \eqref{eq_59} we
conclude that
\begin{eqnarray}
&&\left|\int_0^t \la\tau(u)\nu - S_h(u), 
\frac{\p \gamma^s}{\p t}\ra_{\p\O_j \setminus \G_p}(s) ds\right|\label{eq_60}\\
&&\qquad \le \widehat\epsilon \|\gamma^s(t)\|^2_{1,h} 
+ C\left( \int_0^t \|\gamma^s(s)\|^2_{1,h} ds + h^2
M_1^2\right).\nonumber
\end{eqnarray}
Then add the inequality
\begin{eqnarray*}
\frac{\zeta}{2} \frac{d}{dt} \|\gamma^B\|^2_{0,\O_p} \le \frac{\zeta}{2}
\left(\|\gamma^B\|^2_{0,\O_p} + \|\frac{\p \gamma^B}{\p t}\|^2_{0,\O_p}\right),  
\end{eqnarray*}
to \eqref{eq_48}, integrate the result from $0$ to $T$ 
 and use the bounds \eqref{eq_63}
\eqref{eq_64} \eqref{eq_69} \eqref{eq_62}  and \eqref{eq_60}   to obtain
\begin{eqnarray}
&&\frac12 \left({\mathcal P} \frac{\p \gamma^B}{\p t}, 
\frac{\p \gamma^B}{\p t}\right)_{\O_p}(t)  
+  {\mathcal A}_{\zeta,h}(\gamma^B,\gamma^B)(t)
+ \int_0^t\left( \frac{\eta}{\kappa_0} \frac{\p \gamma^f}{\p t}, \frac{\p
\gamma^f}{\p t} 
\right)_{\O_p}(s) ds \nonumber\\
&&\qquad  + \int_0^t\left\langle {\mathcal D} S_{\G_p}(\frac{\p \gamma^B}{\p t}),
 S_{\G_p}(\frac{\p \gamma^B}{\p t})\right\rangle_{\G_p}(s) ds \label{eq_71}\\
&&\qquad \le \widehat \epsilon \left(\|\gamma^B(t)\|^2_{{\mathcal Z}^h}
+ \int_0^t\left( \frac{\eta}{\kappa_0} \frac{\p \gamma^f}{\p t}, \frac{\p
\gamma^f}{\p t}\right)_{\O_p}(s) ds \right.\nonumber\\
&&\qquad \left.+\int_0^t\left\langle {\mathcal D} S_{\G_p}(\frac{\p \gamma^B}{\p t}),
 S_{\G_p}(\frac{\p \gamma^B}{\p t})\right\rangle_{\G_p}(s) ds  \right) \nonumber\\
&&\qquad + C\left(\int_0^t \left(\|\gamma^B(s)\|_{{\mathcal Z}^h} + 
\|\frac{\p \gamma^B}{\p t}(s)\|_{0,\O_p}\right) ds  + h^2 (M_0^2 +M_1^2)\right)
\end{eqnarray}
Then, use that  ${\mathcal P}$ and ${\mathcal D}$ are positive definite,  
that ${\mathcal A}_{\zeta,h}$ is ${\mathcal Z}^h$-coercive (see \eqref{coercivity_ah}) 
and apply Gronwall's lemma in \eqref{eq_71} to obtain the following estimate:
\begin{eqnarray}
&&\|\frac{\p \gamma^B}{\p t}\|_{L^{\infty}(J,L^2(\O_p)} 
+ \|\gamma^B\|_{L^{\infty}(J,{\mathcal Z}^h}
+ \|\frac{\p \gamma^s}{\p t}\|_{L^{2}(J,L^2(\G_p)} \label{eq_77}\\
&&\hskip2cm  
+ \|\frac{\p \gamma^f\cdot \nu}{\p t}\|_{L^{2}(J,L^2(\G_p)} \le CH^{1/2}( M_0 +
M_1).\nonumber 
\end{eqnarray}
Next, choose $\psi = v^s = v^f =0$, $\varphi = \gamma^H_2$ in \eqref{eq_47}
to get 
\begin{eqnarray}
&&(\text{curl} ~\gamma^E, \gamma^H_2)  + \frac12\frac{d}{dt}(\mu\gamma^H_2,
\gamma^H_2) = \left(\text{curl}~\left[\Pi_h E - E\right],  \gamma^H_2\right)
\label{eq_78}\\
&&\hskip6cm + \left(\mu \frac{\p (P_h H_2 - H_2)}{\p t},  \gamma^H_2\right).\nonumber
\end{eqnarray}
Thus, taking $\psi = \gamma^E$, $\varphi=0$, $V^s = v^f =0$ in \eqref{eq_47}
and using \eqref{eq_78} yields 
\begin{eqnarray}
&&\frac12 \frac{d}{dt}\left[(\varepsilon \gamma^E, \gamma^E) + (\mu  \gamma^H_2,
\gamma^H_2)\right] + (\sigma \gamma^E, \gamma^E) + 
\left\langle \left(\frac{\varepsilon}{\mu}\right)^{1/2}  \gamma^E\cdot \chi, 
\gamma^E\cdot\chi\right\rangle\nonumber\\
&&\qquad  =\left(\frac{\p (\Pi_h E - E)}{\p t},  \gamma^E\right) +
 \left(\text{curl}~\left[\Pi_h E - E\right],  \gamma^H_2\right)\label{eq_80}\\ 
&&\qquad + \left(\mu \frac{\p (P_h H_2 - H_2)}{\p t},  \gamma^H_2\right)
  - \left(L_0
\frac{\eta}{\kappa_0} \frac{\p \gamma^f}{\p t}, \gamma^E\right)_{\O_p} 
+ \left(L_0 \frac{\eta}{\kappa_0} \frac{\p (Q_h u^f- u^f)}{\p t},
\gamma^E\right)_{\O_p}\nonumber\\ 
&&\qquad\qquad+ \left(\sigma\left[\Pi_h E - E\right],  \gamma^E\right) 
- \left(\left[P_h H_2 - H_2\right],  \text{curl}~ \gamma^E\right)\nonumber\\
&&\qquad \qquad 
+ \left\langle \left(\frac{\varepsilon}{\mu}\right)^{1/2}  \left[\Pi_h E - E\right]\cdot \chi, 
\gamma^E\cdot\chi\right\rangle.\nonumber 
\end{eqnarray}
Let us bound each term in the right-hand side of \eqref{eq_80}. First note that 
 note that the assumption that the coefficients $\varepsilon$ and $\mu$ are piecewise
 contants on $\G$, the fact that  $\text{curl}~ \gamma^E\in {\mathcal W}^h$ 
 and the orthogonality property \eqref{orthog_1} that the last two  
 terms in the right-hand
side of  \eqref{eq_80} vanish. 
The other terms can be bounded using the
approximating properties of $\Pi_h$ and $P_h$ in 
\eqref{aprox_pih1}, \eqref{aprox_pih2} and 
\eqref{aprox_ph} as follows.
\begin{eqnarray}
&&\left|\left(\frac{\p (\Pi_h E - E)}{\p t},  \gamma^E\right)\right|
+\left|\left(\text{curl}~\left[\Pi_h E - E\right],  \gamma^H_2\right)\right|
\label{eq_81_82}\\
&&\qquad + \left|\left(\mu \frac{\p (P_h H_2 - H_2)}{\p t},  \gamma^H_2\right)\right|
\nonumber\\
&&\qquad \le C\left(\|\gamma^E_2\|_0^2 +\|\gamma^H_2\|_0^2  
+ h^2 \left( \|\text{curl}~E\|^2_1 + \|\frac{\p E}{\p t}\|^2_1 + 
\|\frac{\p H_2}{\p t}\|^2_1\right) \right).\nonumber. 
\end{eqnarray}
Next use the estimate \eqref{eq_77} to see that 
\begin{eqnarray}
&&\left|\left(L_0 \frac{\eta}{\kappa_0} \frac{\p \gamma^f}{\p t},
\gamma^E\right)\right| \le C\left(\|\gamma^E\|^2_{0,\O_p} 
+  \|\frac{\p \gamma^f}{\p t}\|^2_{0.\O_p}\right)\label{eq_83}\\
&&\hskip3.5cm  \le C\left(\|\gamma^E\|^2_{0,\O_p} +  h (M_0 +
M_1)\right).\nonumber
\end{eqnarray}
Also,
\begin{eqnarray} 
&&\left|\left(L_0 \frac{\eta}{\kappa_0} \frac{\p (Q_h u^f- u^f)}{\p t},
\gamma^E\right)_{\O_p}\right|+ \left|\left(\sigma\left[\Pi_h E - E\right], 
\gamma^E\right)\right|\label{eq_84_85}\\
&&\qquad \le  C\left(\|\gamma^E\|^2_{0} +  h^2 \left(\|E\|^2_1 + 
\|\frac{\p u^f}{\p t}\|^2_{1,\O_p}\right)\right).\nonumber 
\end{eqnarray}
Thus, apply  the bounds \eqref{eq_81_82}, \eqref{eq_83} and \eqref{eq_84_85} in
\eqref{eq_80}  to get the inequality
\begin{eqnarray}
&&\frac12 \frac{d}{dt}\left[(\varepsilon \gamma^E, \gamma^E) + (\mu  \gamma^H_2,
\gamma^H_2)\right] + (\sigma \gamma^E, \gamma^E) \label{eq_86}\\
&&\quad \le  C\left(\|\gamma^E\|^2_{0} + \|\gamma^H_2\|^2_0 +
h^2 \left(\|E\|^2_1+ \|\text{curl}~E\|^2_1 + \|\frac{\p E}{\p t}\|^2_1 +
\|\frac{\p H_2}{\p t}\|^2_1
+ \|\frac{\p u^f}{\p t}\|^2_{1,\O_p}\right)\right. \nonumber\\
&&\qquad \left.+ h( M_0^2 + M_1^2)\right).\nonumber
\end{eqnarray}
Now integrate \eqref{eq_86} from $0$ to $t$ and use \eqref{pos_epsilon_mu} 
  to get 
\begin{eqnarray}
&& \|\gamma^E(t)\|^2_0 + \|\gamma^H_2(t)\|^2 + \int_0^t (\sigma \gamma^E,
\gamma^E)(s) ds\label{eq_87}\\
&&\quad \le  C\left( \|\gamma^E(0)\|^2_{0} + \|\gamma^H_2(0)\|^2_0 + 
\int_0^t \left(\gamma^E(s)\|^2_0 + \|\gamma^H_2(s)\|^2_0\right) ds \right. \nonumber\\
&&\qquad \left.+ C  h^2 (M_1^2 + N_1^2)  +  h( M_0^2 + M_1^2)\right).\nonumber
\end{eqnarray}
Next, choose $E^h(0)$ to be the $\P_h$-projection of $E(0)$ into ${\mathcal V}^h$
and  $H^h_2(0)$ to be the $P_h$-projection of $H^h_2(0)$ into ${\mathcal W}^h$,
so that using the triangle inequality we have  that
\begin{eqnarray}
\|\gamma^E(0)\| \le C h \|E(0)\|_1, \qquad  
\|\gamma^H_2(0)\| \le C h \|\gamma^H_2(0)\|_1.\label{eq_88_1_93}.
\end{eqnarray}
Thus use \eqref{eq_88_1_93} and apply Gronwall's lemma in \eqref{eq_87} to obtain
the estimate
\begin{eqnarray}
\|\gamma^E\|_{L^{\infty}(J,L^2(\O)} 
+ \|\gamma^H_2\|_{L^{2}(J,L^2(\O)}\le C\left( h (N_0 + N_1) + h^{1/2}(M_0 +
M_1)\right). \label{eq_94}
\end{eqnarray}
Finally, using the triangle inequality, the approximating properties 
\eqref{aprox_pih1}, \eqref{aprox_ph} \eqref{approx_rh} and the estimates 
\eqref{eq_77} and \eqref{eq_94} we 
obtain the estimate in \eqref{apriori}. This completes the proof. 
\end{proof}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%  SECTION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



 
\section{The discrete-time finite element procedures for the PSVTM-mode}
Let $L$ a positive integer, $\Delta t = T/L$, $g^n = g(n \Delta t)$. Set 
$\widehat g^{n+1/2} = \frac{g^{n+1} + g^n}{2}$ and  
\begin{eqnarray*}
\p^2 g^n = \frac{g^{n+1} - 2 g^n + g^{n-1}}{\Delta t^2},\ \ 
\p g^n = \frac{g^{n+1} - g^{n-1}}{2 \Delta t},\ \ 
d_t g^n = \frac{g^{n+1} - g^n}{\Delta t}.
\end{eqnarray*}
Our first fully implicit discrete-time procedure combines a Crank-Nicholson scheme for the 
 Maxwell part and a Galerkin  procedure for the Biot-part as follows. 
 Given $\left(E^{h,1}, H^{h,1}_2, u^{s,h,0}, u^{s,h,1}, u^{f,h,0},
 u^{f,h,1}\right)\in  {\mathcal V}^h \times {\mathcal W}^h 
 \times ({\mathcal NC}^h)^4 \times ({\mathcal M}^h)^4$, for $n\ge 1$ compute 
$\left(E^{h,n+1}, H^{h,n+1}_2, u^{s,h,n+1}, u^{s,h,n+1}\right)\in  
{\mathcal Y}^h$ such that
\begin{eqnarray}
&&(\varepsilon ~ d_t E^{h,n}, \psi) +
 (\sigma  \widehat E^{h,n+1/2},\psi) -(\widehat H_2^{h,n+1/2}, \text{curl} ~\psi)\label{eq_103}\\
&& \qquad + \left( L_0 \frac{\eta}{\kappa_0} \p u^{f,h,n}, \psi\right)_{\O_p}
+ \left\langle \left(\frac{\varepsilon}{\mu}\right)^{1/2}  \widehat E^{h,n+1/2}\cdot \chi, 
\psi\cdot\chi\right\rangle=0, \quad \psi\in {\mathcal V}^h,\nonumber \\
&&(\text{curl} \widehat E^{h,n+1/2}, \varphi)  + (\mu d_t H_2^{h,n}, \varphi) = 0, \quad
\varphi\in {\mathcal W}^h,\label{eq_104}\\
&&\left({\mathcal P} \p^2 u^{h,n}, v\right)_{\O_p} 
+ \left( \frac{\eta}{\kappa_0} \p u^{f,h,n}, v^f\right)_{\O_p}
+  {\mathcal A}_h(\frac{u^{h,n+1} + u^{h,n-1}}{2},v) \label{eq_105}\\
&&+ \left\langle {\mathcal D} S_{\G}(\p u^{h,n}),
 S_{\G}(v)\right\rangle_{\G_p} = (F^{n},v), \quad
 v=(v^s,v^f)\in ({\mathcal NC}^h)^2\times {\mathcal M}^h.\nonumber
\end{eqnarray}
To analyze the procedure \eqref{eq_103}-\eqref{eq_105}, choose
$\psi=\widehat E^{h,n+1/2}$ in 
\eqref{eq_103} and $\varphi = \widehat H_2^{h,n+1/2}$ in \eqref{eq_104} and add the
resulting equations to get
\begin{eqnarray}
&&\frac{1}{2 \Delta t}\left[(\varepsilon ~ E^{h,n+1}, E^{h,n+1}) - 
(\varepsilon ~ E^{h,n}, E^{h,n}) + (\mu H_2^{h,n+1}, H_2^{h,n+1}) - 
(\mu H_2^{h,n}, H_2^{h,n})\right]\nonumber\\
&& \qquad + \left( L_0 \frac{\eta}{\kappa_0} \p u^{f,h,n}, E^{h,n+1} \right)_{\O_p}
+ \left\langle \left(\frac{\varepsilon}{\mu}\right)^{1/2}  \widehat E^{h,n+1/2}\cdot \chi, 
E^{h,n+1}\cdot\chi\right\rangle=0.\label{eq_106}
\end{eqnarray} 
Next, choose $v = \p u^{h,n}$ in \eqref{eq_105} and add the resulting equation to
\eqref{eq_106} to obtain 
\begin{eqnarray}
&&\frac{1}{2 \Delta t}\left[(\varepsilon ~ E^{h,n+1}, E^{h,n+1}) - 
(\varepsilon ~ E^{h,n}, E^{h,n}) + (\mu H_2^{h,n+1}, H_2^{h,n+1})\right.\label{eq_110}\\
&&\qquad \left. - (\mu H_2^{h,n}, H_2^{h,n}) 
+ \left({\mathcal P} d_t u^{h,n}, d_t u^{h,n}\right)_{\O_p} - 
  \left({\mathcal P} d_t u^{h,n-1}, d_t u^{h,n-1}\right)_{\O_p}\right.\nonumber\\
&&\qquad \left.  + \frac14 {\mathcal A}_h(u^{h,n+1},  u^{h,n+1}) 
 - \frac14 {\mathcal A}_h(u^{h,n-1},  u^{h,n-1})\right]\nonumber\\
&& \qquad + \Phi(\widehat E^{h,n+1/2}, \p u^{f,h,n}) 
+ (\sigma \widehat E^{h,n+1/2}, \widehat E^{h,n+1/2})_{\O_a}\nonumber\\
&&\qquad + \left\langle \left(\frac{\varepsilon}{\mu}\right)^{1/2}  \widehat E^{h,n+1/2}\cdot \chi, 
E^{h,n+1/2}\cdot\chi\right\rangle + 
\left\langle {\mathcal D} S_{\G}(\p u^{h,n}),
 S_{\G}(\p u^{h,n})\right\rangle_{\G_p} = (F^{n}, \p u^{h,n}),\nonumber
\end{eqnarray} 
where $\Phi(\widehat E^{h,n+1/2}, \p u^{f,h,n})$ is defined in \eqref{def_Phi}. 
Add the inequality 
\begin{eqnarray}
&&\frac{\zeta}{4 \Delta t}\left( \|u^{h,n+1}\|^2_{0.\O_p} -
\|u^{h,n-1}\|^2_{0.\O_p}\right)\nonumber\\
&&\qquad\qquad  \le C \left( \|u^{h,n+1}\|^2_{0.\O_p}
 + \|u^{h,n}\|^2_{0.\O_p} + \|u^{h,n-1}\|^2_{0.\O_p}\right.\label{eq_zeta}\\ 
&&\qquad\qquad \quad \left.  +\|d_t u^{h,n}\|^2_{0.\O_p} + \|d_t
u^{h,n-1}\|^2_{0.\O_p}\right)\nonumber
\end{eqnarray} 
to \eqref{eq_110}, multiply the result by $\Delta t$,  add from $n=1$ to $n=N$
and use \eqref{eq_21} and \eqref{coercivity_ah} to obtain 
\begin{eqnarray}
&&(\varepsilon ~ E^{h,N+1}, E^{h,N+1}) + 
(\mu H_2^{h,N+1}, H_2^{h,N+1})
+ \left({\mathcal P} d_t u^{h,N}, d_t u^{h,N}\right)_{\O_p}\label{eq_113}\\
&&\qquad + \frac{C_2}{4} \left(\|u^{h,N+1}\|^2_{{\mathcal Z}^h} 
+ \|u^{h,N}\|^2_{{\mathcal Z}^h}\right)
+ C_1 \sum_{n=1}^N \left( \|\widehat E^{h,n+1/2}\|^2_{0,\O_p} + \|\p
u^{f,h,n})\|^2_{0,\O_p}\right) \Delta t\nonumber\\
&&\qquad +  \sum_{n=1}^N\left(  (\sigma \widehat E^{h,n+1/2}, \widehat E^{h,n+1/2})_{\O_a} + 
 \left\langle \left(\frac{\varepsilon}{\mu}\right)^{1/2}  \widehat E^{h,n+1/2}\cdot \chi, 
E^{h,n+1}\cdot\chi\right\rangle\right.\nonumber\\
&&\qquad \left.+ \left\langle {\mathcal D} S_{\G}(\p u^{h,n}),
 S_{\G}(\p u^{h,n})\right\rangle_{\G_p}\right)\Delta t\nonumber\\
&&\qquad \le C \left(\|E^{h,1}\|^2_0 +\|H^{h,1}_2\|^2_0 +
\|u^{h,0}\|^2_{{\mathcal Z}^h} + \|u^{h,1}\|^2_{{\mathcal Z}^h}\right.\nonumber\\ 
&&\qquad \left. + \|d_t u^{h,0}\|^2_{0,\O_p} + \sum_{n=1}^N\left(\|F^{n}\|^2_{0,\O_p} + 
+ \|u^{h,n+1}\|^2_{0,\O_p} + \|u^{h,n}\|^2_{0,\O_p} 
+ \|u^{h,n-1}\|^2_{0,\O_p}\right.\right.\nonumber\\ 
&&\qquad \left.\left.+ \|d_t u^{h,n}\|^2_{0,\O_p} + \|d_t u^{h,n-1}\|^2_{0,\O_p} 
 \right) \Delta t\right).\nonumber
\end{eqnarray} 
Next apply Gronwall's lemma in \eqref{eq_113}, use \eqref{pos_epsilon_mu} and
that the matrices ${\mathcal P}$ and ${\mathcal D}$ are positive definite to
conclude that
\begin{eqnarray}
&&\text{max}_{1\le n\le L-1} \left[\|E^{h,n}\|^2_0 + \|H^{h,n}\|^2_0 
+ \|d_t u^{h,n}\|^2_{0,\O_p} 
+ \|u^{h,n}\|^2_{{\mathcal Z}^h}\right]\label{eq_115}\\
&&\qquad + \sum_{n=1}^{L-1} \left( \|\widehat E^{h,n+1/2}\|^2_{0,\O_p} 
+ \|\p u^{f,h,n}\|^2_{0,\O_p} \Delta t\right)\nonumber\\
&&\qquad \sum_{n=1}^{L-1}\left( (\sigma \widehat E^{h,n+1/2}, \widehat E^{h,n+1/2})_{\O_a} + 
 \|\widehat E^{h,n+1/2}\cdot \chi\|^2_{0,\G} \right.\nonumber\\
 &&\qquad\quad  \left. + \|\p u^{s,h,n}\|^2_{0,\G_p} 
 + \|\p u^{f,h,n}\cdot \nu\|^2_{0,\G_p}\right)\Delta t\nonumber\\
&&\qquad \le C \left(\|E^{h,1}\|_0 +\|H^{h,1}_2\|_0 +
\|u^{h,0}\|_{{\mathcal Z}^h} + \|u^{h,1}\|_{{\mathcal Z}^h}
+ \|d_t u^{h,0}\|_{0,\O_p}
\right.\nonumber\\ 
&&\qquad \quad  \left. 
+ \sum_{n=1}^{L-1}
\|F^{n}\|^2_{0,\O_p}\Delta t\right). 
\end{eqnarray}
The estimate \eqref{eq_115} yields uniqueness, existence and unconditional stability for
 the procedure \eqref{eq_103}-\eqref{eq_105}.
  
A second fully implicit in time procedure combines a Backward Euler in time for Maxwell
equations with the implicit procedure used for Biot's equations 
  in \eqref{eq_105}. It can be
formulated as follows: 
Given $\left(E^{h,1}, H^{h,1}_2, u^{s,h,0}, u^{s,h,1}, u^{f,h,0},
 u^{f,h,1}\right)\in  {\mathcal V}^h \times {\mathcal W}^h 
 \times ({\mathcal NC}^h)^4 \times ({\mathcal M}^h)^4$, for $n\ge 1$ compute 
$\left(E^{h,n+1}, H^{h,n+1}_2, u^{s,h,n+1}, u^{s,h,n+1}\right)\in  
{\mathcal Y}^h$ such that
\begin{eqnarray}
&&(\varepsilon ~ d_t E^{h,n}, \psi) +
 (\sigma  \widehat E^{h,n+1},\psi) -(\widehat H_2^{h,n+1/2}, \text{curl}
 ~\psi)\label{eq_153}\\
&& \qquad + \left( L_0 \frac{\eta}{\kappa_0} \p u^{f,h,n}, \psi\right)_{\O_p}
+ \left\langle \left(\frac{\varepsilon}{\mu}\right)^{1/2}  \widehat E^{h,n+1}\cdot \chi, 
\psi\cdot\chi\right\rangle=0, \quad \psi\in {\mathcal V}^h,\nonumber \\
&&(\text{curl} \widehat E^{h,n+1}, \varphi)  + (\mu d_t H_2^{h,n}, \varphi) = 0, \quad
\varphi\in {\mathcal W}^h,\label{eq_154}\\
&&\left({\mathcal P} \p^2 u^{h,n}, v\right)_{\O_p} 
+ \left( \frac{\eta}{\kappa_0} \p u^{f,h,n}, v^f\right)_{\O_p}
+  {\mathcal A}_h(\frac{u^{h,n+1} + u^{h,n-1}}{2},v) \label{eq_155}\\
&&+ \left\langle {\mathcal D} S_{\G}(\p u^{h,n}),
 S_{\G}(v)\right\rangle_{\G_p} = (F^{n},v), \quad
 v=(v^s,v^f)\in ({\mathcal NC}^h)^2\times {\mathcal M}^h.\nonumber
\end{eqnarray}
This procedure can be analyzed with a similar argument to that given for the procedure 
 \eqref{eq_103}-\eqref{eq_105}. The analysis, omitted here for brevity, 
  yields existence, uniqueness and
 unconditional stability for this method. 
 
Our third   discrete-time procedure consists of  Leapfrog scheme for the 
 Maxwell part with the same implicit-time procedure used above 
 for the Biot-part and is formulated  as follows. 
 Given $\left(E^{h,0}, H^{h,1/2}_2, u^{s,h,0}, u^{s,h,1}, u^{f,h,0},
 u^{f,h,1}\right)\in  {\mathcal V}^h \times {\mathcal W}^h 
 \times ({\mathcal NC}^h)^4 \times ({\mathcal M}^h)^4$, for $n\ge 1$ compute 
$\left(E^{h,n}, H^{h,n+1/2}_2, u^{s,h,n+1}, u^{s,h,n+1}\right)\in  
{\mathcal Y}^h$ such that
\begin{eqnarray}
&&(\varepsilon ~ d_t E^{h,n}, \psi) +
 (\sigma  \widehat E^{h,n+1/2},\psi) -(H_2^{h,n+1/2}, \text{curl} ~\psi)\label{eq_116}\\
&& \qquad + \left( L_0 \frac{\eta}{\kappa_0} \p u^{f,h,n}, \psi\right)_{\O_p}
+ \left\langle \left(\frac{\varepsilon}{\mu}\right)^{1/2}  \widehat E^{h,n+1/2}\cdot \chi, 
\psi\cdot\chi\right\rangle=0, \quad \psi\in {\mathcal V}^h,\nonumber \\
&&(\text{curl} \widehat E^{h,n+1/2}, \varphi)  + (\mu \frac{H_2^{h,n+3/2} -
H_2^{h,n+1/2}}{\Delta t}, \varphi) = 0, \quad
\varphi\in {\mathcal W}^h,\label{eq_117}\\
&&\left({\mathcal P} \p^2 u^{h,n}, v\right)_{\O_p} 
+ \left( \frac{\eta}{\kappa_0} \p u^{f,h,n}, v^f\right)_{\O_p}
+  {\mathcal A}_h(\frac{u^{h,n+1} + u^{h,n-1}}{2},v) \label{eq_118}\\
&&+ \left\langle {\mathcal D} S_{\G}(\p u^{h,n}),
 S_{\G}(v)\right\rangle_{\G_p} = (F^{n},v), \quad
 v=(v^s,v^f)\in ({\mathcal NC}^h)^2\times {\mathcal M}^h.\nonumber
\end{eqnarray}
To analyze the procedure \eqref{eq_116}-\eqref{eq_118}, choose $\psi = \widehat
E^{h,n+1/2}$ in \eqref{eq_116} and $\varphi = H_2^{h,n+1/2}$ in \eqref{eq_117}, add the
resulting equations
% and use the identity
%\begin{eqnarray*}
%&&(\mu \frac{H_2^{h,n+3/2} -
%H_2^{h,n+1/2}}{\Delta t}, H_2^{h,n+1/2}) 
%=(\mu \frac{H_2^{h,n+3/2} - H_2^{h,n+1/2}}{\Delta t}, H_2^{h,n+3/2})\\ 
%&&\qquad - (\mu \frac{H_2^{h,n+3/2} - H_2^{h,n+1/2}}{\Delta t}, 
%H_2^{h,n+3/2} - H_2^{h,n+1/2})
%\end{eqnarray*}
to obtain 
%the equation
\begin{eqnarray}
&&\frac{1}{2 \Delta t}\left[(\varepsilon E^{h,n+1}, E^{h,n+1}) - (\varepsilon E^{h,n},
E^{h,n})\right]
+  (\sigma  \widehat E^{h,n+1/2},\widehat E^{h,n+1/2} )\label{eq_120}\\
&&\qquad + (\mu \frac{H_2^{h,n+3/2} - H_2^{h,n+1/2}}{\Delta t},
H_2^{h,n+3/2})\nonumber\\ 
&&\qquad - (\mu \frac{H_2^{h,n+3/2} - H_2^{h,n+1/2}}{\Delta t}, H_2^{h,n+3/2} -
H_2^{h,n+1/2})+ \left( L_0 \frac{\eta}{\kappa_0} \p u^{f,h,n}, 
\widehat E^{h,n+1/2}\right)_{\O_p}\nonumber\\
&&\qquad + \left\langle \left(\frac{\varepsilon}{\mu}\right)^{1/2}  
\widehat E^{h,n+1/2}\cdot \chi, 
E^{h,n+1/2}\cdot\chi\right\rangle=0.\nonumber
\end{eqnarray}
Now add \eqref{eq_120} and the corresponding equation to  \eqref{eq_120} for the case 
$n=n+1$  and use that \eqref{eq_117}
yields
\begin{eqnarray}\label{eq_nueva}
\text{curl} ~\widehat  E^{h,n+1/2} = -\mu \left(\frac{H_2^{h,n+3/2} - H_2^{h,n+1/2}}{\Delta
t}\right)
\end{eqnarray}

\vskip1cm
{\bf Fabio: fijate si las unidades de la ecuacion \eqref{eq_nueva} 
 estan bien, para
obtener esa ecuacion dije lo siguiente: tomo funcion de prueba 

$$
\varphi = \text{curl} ~ \widehat E^{h,n+1/2} + \mu \left(\frac{H_2^{h,n+3/2} - H_2^{h,n+1/2}}{\Delta
t}\right)     
$$

en \eqref{eq_117}  para obtener

$$
\|\text{curl} ~ \widehat E^{h,n+1/2} 
+ \mu \left(\frac{H_2^{h,n+3/2} - H_2^{h,n+1/2}}{\Delta
t}\right)\|^2_0 = 0  
$$

Como la norma es cero, la funcion es cero en $L^2(\O)$ y entonces  

$$
\text{curl} ~\widehat  E^{h,n+1/2} + \mu \left(\frac{H_2^{h,n+3/2} - H_2^{h,n+1/2}}{\Delta
t}\right) = 0 
$$

que es lo que dice arriba en \eqref{eq_nueva}. No les pongo nro para 
no confundirte mas  

Esta es la manera en que aparece este termino (curl, curl)  en el miembro derecho 
de la \eqref{eq_124} debajo, que es el que trae el problema de unidades. }

\vskip1cm 


to get the equation 
\begin{eqnarray}
&&\frac{1}{2 \Delta t}\left[(\varepsilon E^{h,n+2}, E^{h,n+2}) - (\varepsilon E^{h,n},
E^{h,n})\right]
+  (\sigma  \widehat E^{h,n+3/2},\widehat E^{h,n+3/2} ) \label{eq_124}\\
&&\qquad + (\sigma  \widehat E^{h,n+1/2},\widehat E^{h,n+1/2} )
+ (\mu \frac{H_2^{h,n+5/2} - H_2^{h,n+1/2}}{\Delta t}, H_2^{h,n+3/2})\nonumber\\ 
&&\qquad + \left( L_0 \frac{\eta}{\kappa_0} \p u^{f,h,n}, 
\widehat E^{h,n+1/2}\right)_{\O_p}  + \left( L_0 \frac{\eta}{\kappa_0} \p u^{f,h,n+1}, 
\widehat E^{h,n+3/2}\right)_{\O_p}  \nonumber\\
&&\qquad + \left\langle \left(\frac{\varepsilon}{\mu}\right)^{1/2}  
\widehat E^{h,n+3/2}\cdot \chi, E^{h,n+3/2}\cdot\chi\right\rangle
\nonumber\\ 
&&\qquad 
+ \left\langle \left(\frac{\varepsilon}{\mu}\right)^{1/2}  
\widehat E^{h,n+1/2}\cdot \chi, E^{h,n+1/2}\cdot\chi\right\rangle \nonumber\\
&&\qquad \qquad = \left(\text{curl} ~ E^{h,n+1/2}, \frac{\Delta t}{\mu} \text{curl} ~ E^{h,n+1/2}\right).\nonumber
\end{eqnarray}
Next add \eqref{eq_124} to the two equations obtained by choosing $v = \p u^{h,n}$ and 
$v = \p u^{h,n+1}$  in  \eqref{eq_118} to get 
\begin{eqnarray}
&&\frac{1}{2 \Delta t}\left[(\varepsilon ~ E^{h,n+2}, E^{h,n+2}) - 
(\varepsilon ~ E^{h,n}, E^{h,n}) 
+ \left({\mathcal P} d_t u^{h,n+1}, d_t u^{h,n+1}\right)_{\O_p} \right.\label{eq_127}\\
&&\qquad - \left. \left({\mathcal P} d_t u^{h,n-1}, d_t u^{h,n-1}\right)_{\O_p} 
 + \frac14 {\mathcal A}_h(u^{h,n+2},  u^{h,n+2}) 
+ \frac14 {\mathcal A}_h(u^{h,n+1},  u^{h,n+1})\right.\nonumber\\
&&\qquad \left. - \frac14 {\mathcal A}_h(u^{h,n},  u^{h,n})
 - \frac14 {\mathcal A}_h(u^{h,n-1},  u^{h,n-1})\right] 
 +  (\mu \frac{H_2^{h,n+5/2} - H_2^{h,n+1/2}}{\Delta t}, H_2^{h,n+3/2})\nonumber\\
&&\qquad  + \Phi(\widehat E^{h,n+3/2}, \p u^{f,h,n+1})  + 
 + \Phi(\widehat E^{h,n+1/2}, \p u^{f,h,n})
+ (\sigma \widehat E^{h,n+3/2}, \widehat E^{h,n+3/2})_{\O_a}\nonumber\\
&&\qquad + (\sigma \widehat E^{h,n+1/2}, \widehat E^{h,n+1/2})_{\O_a}
+ \left\langle \left(\frac{\varepsilon}{\mu}\right)^{1/2}  \widehat E^{h,n+3/2}\cdot \chi, 
E^{h,n+3/2}\cdot\chi\right\rangle \nonumber\\
&&\qquad + \left\langle \left(\frac{\varepsilon}{\mu}\right)^{1/2}  \widehat E^{h,n+1/2}\cdot \chi, 
E^{h,n+1}\cdot\chi\right\rangle + \left\langle {\mathcal D} S_{\G}(\p u^{h,n+1}),
 S_{\G}(\p u^{h,n+1})\right\rangle_{\G_p}\nonumber\\ 
&&\qquad  + \left\langle {\mathcal D} S_{\G}(\p u^{h,n}),
 S_{\G}(\p u^{h,n})\right\rangle_{\G_p}\nonumber\\ 
&&\qquad  = (F^{n+1}, \p u^{h,n+1}) + (F^{n}, \p u^{h,n}) + 
 \left(\text{curl} ~ E^{h,n+1/2}, \frac{\Delta t}{\mu} \text{curl} ~
 E^{h,n+1/2}\right).\nonumber
\end{eqnarray} 

\vskip1cm
{\bf Fabio: a ver si esto explica lo que pasa y donde meti la pata:

la ecuacion \eqref{eq_127} de arriba esta bien en unidades porque es tomar fuciones de
prueba en las de Maxwell y de Biot y usar 
la \eqref{eq_nueva}, que supongo esta bien y vos vas a chequear. 


Suponiendo eso, pasamos ahora a la desigualdad 
\eqref{eq_131} que sigue , donde fijate por ejemplo que 
la constante $C_1$ tiene unidades,  basta mirar la definicion de la $C_1$, de lo contrario
no se podrian sumar  los  terminos 

$$\|\widehat E^{h,n+1/2}\|^2_{0,\O_P}$$

 y 

$$
\|\p u^{f,h,n}\|^2_{0,\O_p}
$$

Por lo tanto lo mismo pasa me parece en el mienbro derecho de la \eqref{eq_131}, para
poder sumar por ej.  

$$\|E^{h,2}\|^2_0$$

con
 
$$(\mu H_2^{h,3/2}, H_2^{h,1/2})$$

o dentro de la sumatoria , sumar 

$$\|F^n\|^2_{0,\O_p}$$  

con 

$$\|u^{h,n+1}\|^2_{0,\O_p}$$

la constante (que ahora le pongo nombre $C_8$ debe tener unidades, las que correspondan
para poder sumar todos esos terminos antes del ultimo termino, 
que esta sin tocarse todavia. 
 Por lo tanto, para poder incluir ese ultimo termino dentro de la sumatoria 
 debemos pedir que 

$$ 
\Delta t \dfrac{C_4^2 h^{-2} }{2 \varepsilon^* \mu_*}\le C_8  
$$

y deberia pasar que $C_8$ tenga las unidades que hacen falta.

Por favor si te parece esto bien, hace el laburo de estimar la constante $C_8$ ya que la
necesitamos, quedara seguro en funcion de las cotas de los coeficientes de las ecuaciones,
como se hico con la $C_1$.

Yo arreglo esto debajo suponiendo que lo que dije esta bien.
}
\vskip1cm




 
Next, add  \eqref{eq_zeta} and the corresponding inequality for $n=n+1$  
 to \eqref{eq_127}, multiply the resulting inequality by $\Delta t$ and sum from $n=1$
 to $n=N$ and use \eqref{eq_21} and \eqref{coercivity_ah}  to get
\begin{eqnarray}
&&\frac12\left[(\varepsilon E^{h,N+2}, E^{h,N+2}) +  (\varepsilon E^{h,N+1},
E^{h,N+1})\right]  
+ \left({\mathcal P} d_t u^{h,N+1}, d_t u^{h,N+1}\right)_{\O_p} \nonumber\\
&&\quad + \left({\mathcal P} d_t u^{h,N}, d_t u^{h,N}\right)_{\O_p}
 + \frac{C_2}{4} \left( \|u^{h,N+2}\|^2_{{\mathcal Z}^h} 
+ \|u^{h,N+1}\|^2_{{\mathcal Z}^h} + \|u^{h,N}\|^2_{{\mathcal Z}^h}\right)\nonumber\\
&&\quad  +  \sum_{n=1}^{N+1} \left[ C_1 \left(\|\widehat E^{h,n+1/2}\|^2_{0,\O_P} 
 + \|\p u^{f,h,n}\|^2_{0,\O_p} \right) + \|\widehat E^{h,n+1/2}\|_{0,\O_a}\right] \Delta t \label{eq_131}\\
&&\quad + \sum_{n=1}^{N+1}\left(\left\langle \left(\frac{\varepsilon}{\mu}\right)^{1/2}  
\widehat E^{h,n+1/2}\cdot \chi, \widehat E^{h,n+1/2}\cdot\chi\right\rangle
\right.\nonumber\\
&&\qquad \qquad \left. +  \left\langle {\mathcal D} S_{\G}(\p u^{h,n}),
 S_{\G}(\p u^{h,n})\right\rangle_{\G_p}\right) \Delta t
% \nonumber\\
%&&\quad  
+  (\mu H_2^{h,N+5/2}, H_2^{h,N+3/2})\nonumber\\
&&\qquad\qquad  \le C_8 \left[\|E^{h,2}\|^2_0 + \|E^{h,1}\|^2_0 + (\mu H_2^{h,3/2}, H_2^{h,1/2}) 
+ \|d_t u^{h,0}\|^2_{0,\O_p} \right.\nonumber\\
&&\qquad\qquad \left.+ \|d_t u^{h,1}\|^2_{0,\O_p} 
+ \|u^{h,0}\|^2_{{\mathcal Z}^h} +  \|u^{h,1}\|^2_{{\mathcal Z}^h} 
+ \|u^{h,2}\|^2_{{\mathcal Z}^h}\right.\nonumber\\
&&\qquad\qquad \left.+ \sum_{n=1}^N \left(\|F^n\|^2_{0,\O_p} + \|F^{n+1}\|^2_{0,\O_p} 
+ \left({\mathcal P} d_t u^{h,n+1}, d_t u^{h,n+1}\right)_{\O_p}\right. \right.\nonumber\\
&&\qquad\qquad \left. \left.+ \left({\mathcal P} d_t u^{h,n}, d_t u^{h,n}\right)_{\O_p}
+ \left({\mathcal P} d_t u^{h,n-1}, d_t u^{h,n-1}\right)_{\O_p}\right. \right.\nonumber\\  
&&\qquad\qquad \left. \left.+ \|u^{h,n+1}\|^2_{0,\O_p}+ \|u^{h,n}\|^2_{0,\O_p}\right)
\Delta t\right]
\nonumber\\
&&\qquad\qquad +\sum_{n=1}^N  \Delta t\left( \frac{1}{\mu}\text{curl} ~ E^{h,n+1/2}, \text{curl} ~
 E^{h,n+1/2}\right) \Delta t.\nonumber
\end{eqnarray}

\vskip1cm
{\bf Fabio: lo que dije en el mail, mas claro aqui: en la desigualdad de arriba, las
unidades de 

$$\frac12\left[(\varepsilon E^{h,N+2}, E^{h,N+2}) +  (\varepsilon E^{h,N+1},
E^{h,N+1})\right]
$$ 
deberian ser las mismas que las de }

$$
\sum_{n=1}^N  \Delta t\left( \frac{1}{\mu}\text{curl} ~ E^{h,n+1/2}, \text{curl} ~
 E^{h,n+1/2}\right) \Delta t.
$$
\vskip1cm

Let us bound the last term in the right-hand side of  \eqref{eq_131}. By the
quasiregularity assumption on the finite element partition ${\mathcal T}^h$, the exists
a positive constant $C_4$ ($C_4 = \sqrt 48$) such that
\begin{eqnarray}
\hskip3cm \|\text{curl} ~ E^{h,n+1/2}\|_0 \le C_4 h^{-1} \|E^{h,n+1/2}\|_0.\label{eq_133}
\end{eqnarray}
Thus, using \eqref{pos_epsilon_mu},
\begin{eqnarray}
&&\left(\frac{1}{\mu}\text{curl} ~ E^{h,n+1/2}, \text{curl} ~
E^{h,n+1/2}\right)\label{eq_133-1}\\
&&\hskip2cm  \le \dfrac{C_4^2 h^{-2} }{2 \varepsilon^* \mu_*}
 \left( \left(\varepsilon E^{h,n+1}, E^{h,n+1}\right) +\left(\varepsilon E^{h,n+1},
 E^{h,n+1}\right) \right).\nonumber
\end{eqnarray}
Thus we conclude that  
\begin{eqnarray}
&&\sum_{n=1}^N  \Delta t\left( \frac{1}{\mu}\text{curl} ~ E^{h,n+1/2}, \text{curl} ~
 E^{h,n+1/2}\right) \Delta t \label{eq_135-11}\\
&&\qquad \quad \le \sum_{n=1}^N  \Delta t \dfrac{C_4^2 h^{-2} }{2 \varepsilon^* \mu_*}
\left( \left(\varepsilon E^{h,n+1}, E^{h,n+1}\right) +\left(\varepsilon E^{h,n+1},
 E^{h,n+1}\right) \right)\Delta t \nonumber\\
 &&\qquad \quad \le  \sum_{n=1}^N  \left((\varepsilon E^{h,n+1}, E^{h,n+1}) 
+ (\varepsilon E^{h,n}, E^{h,n})\right)\Delta t,\nonumber
\end{eqnarray}
provided 
\begin{eqnarray}\label{stabil_0}
\Delta t \dfrac{C_4^2 h^{-2} }{2 \varepsilon^* \mu_*}\le C_8
\end{eqnarray}
or equivalently,  
assume  the stability constraint
\begin{eqnarray}
\hskip4cm \Delta t \le  C_8 \dfrac{2 \varepsilon^* \mu_* h^2}{C_4^2}\label{stabil_1}.
\end{eqnarray}

\vskip1cm
{\bf Fabio: la desigualdad \eqref{stabil_1} es de la forma clasica  
$\Delta t \le C h^2$ valida para metods explicitos. ES posible que 
las unidades en la \eqref{eq_135-11} esten bien (yo creo que si, al menos 
uno comienza con la  \eqref{eq_133-1} que tiene bien las unidades). }

\vskip1cm

Thus we conclude that  
\begin{eqnarray}
&&\sum_{n=1}^N  \Delta t\left( \frac{1}{\mu}\text{curl} ~ E^{h,n+1/2}, \text{curl} ~
 E^{h,n+1/2}\right) \Delta t \label{eq_135-1}\\
&&\qquad \quad \le  C_8 \sum_{n=1}^N  \left((\varepsilon E^{h,n+1}, E^{h,n+1}) 
+ (\varepsilon E^{h,n}, E^{h,n})\right) \Delta t.\nonumber
\end{eqnarray}
Next we analyze the  last term in the left-hand side of \eqref{eq_131}. 
First, right
this term in the form
\begin{eqnarray}
&&(\mu H_2^{h,N+5/2}, H_2^{h,N+3/2}) = 
(\mu \left[H_2^{h,N+5/2}- H_2^{h,N+3/2})\right], H_2^{h,N+3/2}) \label{eq_135-2}\\
&&\hskip4cm + (\mu H_2^{h,N+3/2}, H_2^{h,N+3/2})  \equiv T_1 + T_2.\nonumber
\end{eqnarray}
Note that from \eqref{eq_117} and  \eqref{eq_133} 

\vskip1cm
{\bf Fabio aca uso  lo de siempre, que para el $\epsilon$ que yo quiera, } 

$a b \le \frac12 (\epsilon a^2 + \frac{1}{\epsilon} b^2)$

\vskip1cm

\begin{eqnarray}
&&|T_1| = \left| - \Delta t \left( \text{curl}~ E^{h,N+3/2},
H_2^{h,N+3/2}\right)\right|\label{eq_135-3}\le \Delta t  \|\text{curl}~ E^{h,N+3/2}\|_0
\|H_2^{h,N+3/2}\|_0\nonumber \\
&&\quad \le \dfrac{\mu_*}{2} \|H_2^{h,N+3/2}\|^2_0 
+ (\Delta  t)^2 \dfrac{C_4^2 h^{-2} }{2\mu_*} \|E^{h,N+3/2}\|^2_0\nonumber\\
&&\quad \le \dfrac{\mu_*}{2} \|H_2^{h,N+3/2}\|^2_0  
+ (\Delta  t)^2 \dfrac{C_4^2 h^{-2} }{4 \varepsilon^* \mu_*}
\left( \left(\varepsilon E^{h,n+1}, E^{h,n+1}\right) +\left(\varepsilon E^{h,n+1},
 E^{h,n+1}\right)\right).\nonumber
\end{eqnarray}
Assume that 
\[
(\Delta  t)^2 \dfrac{C_4^2 h^{-2} }{4 \varepsilon^* \mu_*} \le \frac14
\]
so that, in addition to \eqref{stabil_1}, $\Delta t$ satifies the constraint
\begin{eqnarray}
\hskip4cm \Delta t \le \dfrac{(\mu_* \varepsilon^*)^{1/2}}{C_4} h \label{stabil_2}.
\end{eqnarray}

\vskip1cm
{\bf Fabio: esto parece estar bien, habia un delta t dividiendo que estaba mal


Lo que estoy haciendo es hacer que este termino sea la mitad del termino

$$
\frac12\left[(\varepsilon E^{h,N+2}, E^{h,N+2}) +  (\varepsilon E^{h,N+1},
E^{h,N+1})\right]
$$

por eso pido el $\frac14$.  En otras palabras, lo que estoy usando es que 
$(\Delta  t)^2 \dfrac{C_4^2 h^{-2} }{4 \varepsilon^* \mu_*}$ es adimensional, 
si no  fuera asi no lo puedo restar al $\frac12$.

Si esto esta bien, veamos como arreglar  o entender la primera condicion de
estabilidad  \eqref{stabil_1}, que es practicamente el cuadrado de esta salvo un
factor de 2. Donde carajo esta el error en la derivacion de la  \eqref{stabil_1}? }
\vskip1cm 

Then from \eqref{eq_135-2} and \eqref{eq_135-3}, since 
\[
|T_2| \ge \mu_*\|H_2^{h,N+3/2}\|^2_0,
\]
 
\begin{eqnarray}
&&(\mu H_2^{h,N+5/2}, H_2^{h,N+3/2}) \ge  \dfrac{\mu_*}{2} \|H_2^{h,N+3/2}\|^2_0
\label{eq_136-1}\\
&&\hskip3cm  - \frac14 \left[ \left((\varepsilon E^{h,N+2}, E^{h,N+2}) 
+ (\varepsilon E^{h,N+1}, E^{h,n+1})\right)\right].\nonumber
\end{eqnarray}
Using \eqref{eq_135-1} and \eqref{eq_136-1} in \eqref{eq_131} we obtain
\begin{eqnarray}
&&\frac14\left[(\varepsilon E^{h,N+2}, E^{h,N+2}) +  (\varepsilon E^{h,N+1},
E^{h,N+1})\right] +  \dfrac{\mu_*}{2} \|H_2^{h,N+3/2}\|^2_0 \label{eq_137}\\
&&\quad + \left({\mathcal P} d_t u^{h,N+1}, d_t u^{h,N+1}\right)_{\O_p} 
 + \left({\mathcal P} d_t u^{h,N}, d_t u^{h,N}\right)_{\O_p}\nonumber\\
&&\quad  + \frac{C_2}{4} \left( \|u^{h,N+2}\|^2_{{\mathcal Z}^h} 
+ \|u^{h,N+1}\|^2_{{\mathcal Z}^h} + \|u^{h,N}\|^2_{{\mathcal Z}^h}\right)\nonumber\\
&&\quad  +  \sum_{n=1}^{N+1} \left[ C_1 \left(\widehat E^{h,n+1/2}\|^2_{0,\O_P} 
 + \|\p u^{f,h,n}\|_{0,\O_p} \right) + \|\widehat E^{h,n+1/2}\|_{0,\O_a}\right] \Delta t \\
&&\quad + \sum_{n=1}^{N+1} \left(\left\langle \left(\frac{\varepsilon}{\mu}\right)^{1/2} 
 \widehat E^{h,n+1/2}\cdot \chi, \widehat E^{h,n+1/2}\cdot\chi\right\rangle
 \right.\nonumber\\
&&\qquad\qquad  \left.  +  \left\langle {\mathcal D} S_{\G}(\p u^{h,n}),
 S_{\G}(\p u^{h,n})\right\rangle_{\G_p}\right) \Delta t\nonumber\\
&&\qquad \le C\left[\|E^{h,2}\|^2_0 + \|E^{h,1}\|^2_0 + \|H_2^{h,1/2}\|^2_0 +
\|H_2^{h,3/2}\|^2_0  
+ \|d_t u^{h,0}\|^2_{0,\O_p} \right.\nonumber\\
&&\qquad \left.+ \|d_t u^{h,1}\|^2_{0,\O_p} 
+ \|u^{h,0}\|^2_{{\mathcal Z}^h} +  \|u^{h,1}\|^2_{{\mathcal Z}^h} 
+ \|u^{h,2}\|^2_{{\mathcal Z}^h}\right.\nonumber\\
&&\qquad \left. \sum_{n=1}^N \left(\|F^n\|^2_{0,\O_p} + \|F^{n+1}\|^2_{0,\O_p} 
+  (\varepsilon E^{h,n+2}, E^{h,n+2})\right.\right.\nonumber\\
&&\qquad \quad  \left.\left. + (\varepsilon E^{h,n+1}, E^{h,n+1}) 
+ \left({\mathcal P} d_t u^{h,n+1}, d_t u^{h,n+1}\right)_{\O_p}\right. \right.\nonumber\\ 
&&\qquad\quad \left. \left. + \left({\mathcal P} d_t u^{h,n}, d_t u^{h,n}\right)_{\O_p}
 +  \left({\mathcal P} d_t u^{h,n-1}, d_t u^{h,n-1}\right)_{\O_p} \right. \right.\nonumber\\  
&&\qquad\quad \left. \left.+ \|u^{h,n+1}\|^2_{0,\O_p}+ \|u^{h,n}\|^2_{0,\O_p}\right)
 \Delta t\right].\nonumber
\end{eqnarray}
Next apply the discrete Gronwall's lemma in \eqref{eq_127},  
use  the positive definitess of the matrices  ${\mathcal P}$ and 
${\mathcal D}$ and  \eqref{pos_epsilon_mu}  
to derive the estimate
\begin{eqnarray}
&&\text{max}_{1\le n\le L-1} \left[\|E^{h,n}\|^2_0 + \|H_2^{h,n+1/2}\|^2_0 
+ \|d_t u^{h,n}\|^2_{0,\O_p} 
+ \|u^{h,n}\|^2_{{\mathcal Z}^h}\right]\label{eq_138}\\
&&\qquad + \sum_{n=1}^{L-1} \left( \|\widehat E^{h,n}\|^2_{0,\O_p} 
+ \|\p u^{f,h,n}\|^2_{0,\O_p} \Delta t\right)\nonumber\\
&&\qquad \sum_{n=1}^{L-1}\left( (\sigma \widehat E^{h,n+1/2}, \widehat E^{h,n+1/2})_{\O_a} + 
 \|\widehat E^{h,n+1/2}\cdot \chi\|^2_{0,\G} \right.\nonumber\\
 &&\qquad\quad  \left. + \|\p u^{s,h,n}\|^2_{0,\G_p} 
 + \|\p u^{f,h,n}\cdot \nu\|^2_{0,\G_p}\right)\Delta t\nonumber\\
&&\qquad \le C \left(\|E^{h,1}\|_0 + \|E^{h,2}\|_0 +\|H^{h,1/2}_2\|_0
+\|H^{h,3/2}_2\|_0 +
\|u^{h,0}\|_{{\mathcal Z}^h} + \|u^{h,1}\|_{{\mathcal Z}^h} \right.\nonumber\\
&&\qquad \quad \left. + \|u^{h,2}\|_{{\mathcal Z}^h} 
+ \|d_t u^{h,0}\|_{0,\O_p} + \|d_t u^{h,1}\|_{0,\O_p}
+ \sum_{n=1}^{L-1} \|F^{n}\|^2_{0,\O_p}\Delta t\right).\nonumber 
\end{eqnarray}
As a consequence of the estimate \eqref{eq_138} we have uniqueness and existence for
the procedure \eqref{eq_116}-\eqref{eq_118}, provided the time step $\Delta t$
satisfies the stability conditions \eqref{stabil_1} and \eqref{stabil_2}. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%  SECTION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



 
\section{The triangular case for the PSVTM-mode}
Let ${\mathcal T}^h$ be a quasiregular partition of $\O$ into triangles $\O_j$ of
diameter bounded by $h$. Set
\begin{eqnarray*}
R_1 = \{U : u = (a + b~ x_2, c + d ~x_1), a,b,c \text{constants}\},
\end{eqnarray*}
and let us define the spaces to approximate the electric and magnetic 
fields  ${\mathcal V}^h$ and  ${\mathcal V}^h$ as  
\begin{eqnarray*}
&&{\mathcal V}^h = \{\psi\in H(\text{curl},\O): \psi|_{\O_j}\in R_1(\O_j) \ \forall
j\},\\
&&{\mathcal W}^h = \{\varphi\in L^2(\O): \varphi|_{\O_j}\in P_0(\O_j) \ \forall j\}.
\end{eqnarray*}
The local degrees of fredom for ${\mathcal V}^h$  and ${\mathcal W}^h$
 can be taken to be the values of the
tangential components at the mid-points of each edge of the triangle $\O_j$ and 
 the values at the centroids of
$\O_j$, respectively.
 
Next,  if
\[
{\mathcal NC}^h_j = P_1(\O)_j),
\]
the space ${\mathcal NC}^h$ to approximate each component of the solid displacement in$\O_p$ is
defined as in the rectangular case as 
\begin{eqnarray*}
&&{\mathcal NC}^{h} = \{ v \, :  v_j=  v|_{\O_j}\,  \in 
{\mathcal NC}^h_j, ~ 
v_j(\xi_{jk}) =v_k(\xi_{jk})\,\,  \forall (j,k)\}, 
\end{eqnarray*}
with the local degrees of freedom being the values at the mid-point of the edges of
each $\O_j$.

Finally if 
\begin{eqnarray*}
S_1 = \{U : u = (a + b~ x_1, c + d ~x_2), a,b,c \text{constants}\},
\end{eqnarray*}
the space to approximate the fluid displacement vector is 
\begin{eqnarray*}
&&{\mathcal M}^h = \{v^f\in H(\text{div},\O_p): v^f|_{\O_j}\in S_1(\O_j) \ \forall
j\}.
\end{eqnarray*}
The local degrees of fredom for ${\mathcal M}^h$ are the 
values of normal  components at the mid-points of each edge of the triangle $\O_j$.

The projection operators $\Pi_h, P_h, Q_h, R_h$ and $S_h$ are defined as in the
rectangular case, with identical properties.

Next, setting 
\[
{\mathcal Y}^h = {\mathcal V}^h\times {\mathcal W}^h \times ({\mathcal NC}^h)^2 \times
{\mathcal M}^h
\] 
the definition of the  continuous and discrete-time finite element procedures remain
unchanged, with identical  {\it apriori} error estimates and stability results.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%  SECTION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{3D seismoelectric modeling. 3D rectangular elements.}
Consider  a 3D-rectangular  domain $\O = \O_a \cup \O_p$ 
 where $\O_a$ and $\O_p$ 
are  associated with the air and subsurface poroelastic  (disjoint)  parts 
of $\O$, respectively. Consider the solution of 
\begin{eqnarray} 
&& \varepsilon \frac{\p E}{\p t}  +  \sigma E - \nabla\times  H + L_0 \frac{\eta}{\kappa_0} \frac{\p
u^f}{\p t}  = 0,\quad \O,\label{eq1}\\
&&\nabla\times   E +   \frac{\p  H}{\p t} = 0, \quad \O, 
\label{eq2}\\
&&  \rho_b \frac{\p^2 u^s}{\p t^2}  + \rho_f \frac{\p^2 u^f}{\p t^2}  
 -  \nabla\cdot \tau(u) = F^{(s)},\quad \O_p,
\label{eq3}\\
&& \rho_f \frac{\p^2 u^s}{\p t^2} +  m  \frac{\p^2 u^f}{\p t^2}  + 
  \frac{\eta}{\kappa_0} \ \frac{\p u^f}{\p t} u^f  
 + \nabla p_f = F^{(f)},\quad \O_p,\label{eq4}
\end{eqnarray}
where $\tau(u), p_f(u)$ are given by \eqref{mod.1e} and \eqref{mod.1f}, respectively.
We will solve \eqref{eq1}-\eqref{eq4} with the boundary conditions
\begin{eqnarray}
&&  \varepsilon^{1/2} P_{\chi} E  + \mu^{1/2} \nu\times H = 0, \quad \text{on} \quad 
\G,\label{eq8}\\
&& -{\mathcal G}_{\G_p}(u)  =  {\mathcal D}S_{\G_p}(\frac{\p u}{\p t}), 
\quad \text{on}\quad \G_p,\label{eq10}
\end{eqnarray}
the free surface condition
\begin{eqnarray}
&& -{\mathcal G}_{\G_p}(u)  = 0, 
\quad \text{on}\quad \G_{a,p},\label{eq8a}
\end{eqnarray}
and the initial conditions
\begin{eqnarray}
&&E(t=0) = E_0, \quad H(t=0) = H_{0} 
\label{icond_3d}\\
&& u^s(t=0) = u^s_0, \quad u^f(t=0) = u^f_0,\nonumber\\
&& \frac{\p u^s}{\p t}(t=0) = u^s_1, \quad \frac{\p u^f}{\p
t}(,t=0) = u^f_1.
\nonumber
\end{eqnarray}
Here 
\[
P_{\chi} E = E - \nu (\nu \cdot E) = -\nu\times(\nu\times E)
\]
is the 3-D orthogonal projection of the trace of E into the tangential plane
perpendicular to the normal vector $\nu$ and
\begin{subeqnarray*}
{\mathcal G}_{\G_s}(u) &=& \bigg(\tau(u) \nu \cdot \nu, \tau(u)\nu\cdot \chi^1,
\tau(u)\nu\cdot \chi^2,  p_f(u)\bigg)^t,\\
S_{\G_s}(u) &=& \left(u^s\cdot\nu,u^s\cdot\chi^1,u^s\cdot\chi^1, u^f\cdot\nu\right)^t,
\end{subeqnarray*} 
where , where $^t$ denotes the transpose,  $\nu$ is  the unit outer normal 
on $\G_s$  and 
$\chi^1, \chi^2$ are two   unit tangents on $\G_s$ such that 
$\{\nu, \chi^1, \chi^2\}$ is an orthonormal set on $\G_p$.
The matrices $4\times 4$ matrices  ${\mathcal R}, {\mathcal M}$ in the definition of the matrix 
${\mathcal D}$ in \eqref{eq10} are defined as the obvious extension to the 3D case of
the defining equations in the PSVTM case (simply add a row and a 
column containing all
zeroes and either $b$ or $G$ in the main diagonal as needed).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%  SECTION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{A weak formulation}
First recall the integration by parts formula \cite{dsheenb}
\begin{eqnarray}\label{by_parts3}
&&(\nabla\times U,  V) - (U, \nabla\times V)=\la\nu\times U, V\ra,\,\quad
\forall~ U, V \in H({\text{curl}},\O).
\end{eqnarray}
Thus, multiply \eqref{eq1} by $\psi\in H({\text{curl}},\O)$, integrate in $\O$ 
 and use \eqref{by_parts3}. Also multiply \eqref{eq2} by $\varphi\in [L^2(\O)]^3$ 
 and integrate. Finally, multiply  \eqref{eq3} by $v^s\in [H^1(\O)]^3$ and 
 \eqref{eq4} by $v^f\in H({\text{div}},\O)$, add the resulting equations, 
  integrate in $\O_p$ and use \eqref{eq10}. WE obtain the weak formulation: find $(E, H
  , u^s, u^f)\in  H({\text{curl}},\O)\times [L^2(\O)]^3 
  \times [H^1(\O)]^3\times H({\text{div}},\O)$ such that 
\begin{eqnarray}
&& (\varepsilon \frac{ \p E}{\p t}, \psi)
 + (\sigma E, \psi) -(H, \nabla\times  \psi)  
 + \left( L_0 \frac{\eta}{\kappa_0} \frac{\p  u^f}{\p t}, \psi\right)_{\O_p}
 \label{eq14}\\
&&\hskip3cm + \left\langle \left( \frac{\varepsilon}{\mu}\right)^{1/2} P_{\chi} E,
P_{\chi}\psi \right\rangle
 = 0, \quad \psi \in H(\text{curl},\O), 
\nonumber \\
&&\qquad (\nabla\times E, \varphi) + (\mu \frac{\p H}{\p t}, \varphi) = 0, \quad \varphi
\in [L^2(\O)]^3.\label{eq15}\\
&& ({\mathcal P} \frac{\p^2 u}{\p t^2}, v)_{\O_p} 
+ ( \frac{\eta}{\kappa_0} \frac{\p u^f}{\p t}, v^f)_{\O_p} 
+ {\mathcal A}(u,v) + \left\langle {\mathcal D} S_{\G_p}(\frac{\p u}{\p t}),
S_{\G_p}(v)\right\rangle_{\G_p}\label{eq16}\\
&&\hskip1cm   = (F,v)_{\O_p} , \quad v=(v^s, v^f)\in 
   [H^1(\O_p)]^3 \times H(\text{div},\O_p)\nonumber. 
\end{eqnarray}
In \eqref{eq16} the matrix  ${\mathcal P}$ is defined as in \eqref{def_P} and 
 ${\mathcal A}(u,v)$ is defined as in \eqref{def_A} provided the symmetric $7\times 7$
 matrix  ${\mathcal M} = (m_{ij})$
 is defined with entries $m_{11} = m_{22}= m_{33} = \lambda_c + 2G, 
 m_{12} = m_{13} = \lambda_c, m_{14}= \alpha K_{av}, m_{15} = m_{16} = m_{17} =0, 
 m_{23} = \lambda_c, m_{24}= \alpha K_{av}, m_{25} = m_{26} = m_{27} =0, 
 m{34} = \alpha K_{av}, m_{35} = m_{36} = m_{37} =0, m_{44} =  K_{av}, 
 m_{45} = m_{46} = m_{47} =0, m_{55}= m_{66} = m_{77} = G, 
 m_{56}=m_{57} = m_{67}= 0$. Also, define 

$\tilde\varepsilon(u) = \left(
\e_{11}(u^{(s)}), \e_{22}(u^{(s)}), \e_{33}(u^{(s)}), \nabla\cdot u^{(f)}, 
\e_{12}(u^{(s)},\e_{13}(u^{(s)}, \e_{23}(u^{(s)}\right)^t$.
  
Now equations \eqref{eq14}-\eqref{eq16} are formally identical to equations 
\eqref{eq_13}-\eqref{eq_15} for the 2D case and the argument for uniqueness and
existence  of the
continuos problem  can be repeated to demonstrate that the statement in 
theorem  \thmref{Theorem1}
remains valid in the 3D case.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%  SECTION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{A continuous-time finite element procedure}



 \begin{thebibliography}{}

\bibitem{pride96} S. R. Pride and M. W. Haartsen,
``Electroseismic wave properties," {  J. Acoust. Soc. Amer.}, { 100 (3)},  
1301 (1996).

\bibitem{mikha97} O. V. Mikhailov, M. W. Haartsen and M.  N. Toksoz,
``Electroseismic investigation of the sallow subsurface: Field measurements and
numerical modeling," { Geophysics}, { 62}, 97 (1997).

\bibitem{mikha2000} O. V. Mikhailov, J. Queen and M.  N. Toksoz,
``Using borehole electroseismic measurements to detect and characterize fractured (permeable) zones,"  {Geophysics}, { 65}, 1098 (2000).

\bibitem{thompson05a} A. H. Thompson, 
``Electromagnetic-to-seismic conversion: Successful developments suggest viable applications in exploration and production," 75th SEG Annual Meeting, 
Expanded Abstracts, Houston,  (2005).

\bibitem{hornbostel07} S. Hornbostel and  A. H. Thompson,
``Waveform design for electroseismic exploration," {  Geophysics}, { 72 (2)}, Q1 (2007)

\bibitem{thompson93} A. H. Thompson and G. Gist,
``Geophysical applications of electrokinetic conversion," {  The Leading Edge}, { 12}, 
1169 (1993).


\bibitem{pride_94}
S. R. Pride,  ``Governing equations for the coupled  electromagnetics and acoustics  of porous media," { Physics Review B}, { 50}, 15678 (1994),

\bibitem{block06} G. I. Block and J. G. Harris, 
``Conductivity dependence of seismoelectric wave phenomena in 
fluid-saturated sediments," { J. Geophys. Res.}, 
{ 111}, B01304 (2006), DOI:10.1029/2005JB003798.

\bibitem{haines_pride_06} S. H. Haines  and S. R. Pride, 
``Seismoelectric numerical modeling on a grid,"  
{Geophysics}, 
{ 71 (6)},  N57 (2006)

\bibitem{biot56b} M. A. Biot,
``Theory of propagation of elastic waves in a fluid-saturated 
porous solid. I. Low frequency range," { J. Acoust. Soc. Amer.}, 
{ 28},  168 (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. Amer.}, {28},  179 (1956).

\bibitem{haartsen_pride_97} M. H. Haartsen and S. R. Pride, 
``Electroseismic waves from point sources in  layered media,"  
{ J. Geophys. Res.}, 
{102}, (NO. B11),  24745 (1997).

\bibitem{han_01} Q. Han and  Z. Wang,
``Time-domain simulation of SH-wave-induced electromagnetic field in
heterogeneous porous media: A fast finite element algorithm,"{ Geophysics}, 
{ 66 (2)},  448 (2001)

\bibitem{pain2005} C. C. Pain, J. H. Saunders, M. H. Wortington, J. M. Singer, 
W. Stuart-Bruges, G. Mason and A. Goddard,
``A mixed finite element method for solving the poroelastic Biot equations with
electrokinetic coupling,"  { Geophys. J. Int.}, { 160}, 592 (2005).


\bibitem{bwhite_05}
B. S. White,
``Asymptotic theory of electroseismic prospecting,"
{ SIAM J. Appl. Math.}, { 65} (4), 1443 (2005).

\bibitem{bwhite_06}
B. S. White and M. Zhou,
``Electroseismic prospecting in layered media,"
 { SIAM J. Appl. Math.}, { 67} (1), 69 (2006).

\bibitem{rath} P. A. Raviart   and J. M. Thomas, 
``A mixed finite element
method for second order elliptic problems," in  Mathematical Aspects of
the Finite Element Method, Lecture Notes in Mathematics, volume 606,
 { Springer-Verlag, Berlin, New York}, I. Galligani and E.
Magenes, eds.,  (1977), p.292.

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


\bibitem{pmonk_93} P. B. Monk,  
``An analysis of Nedelc method for the spatial discretization of Maxwell 
equations" 
{ J. Comput. Appl. Math}, 
{ 47},   101 (1993).

\bibitem{pmonk_94} P. B. Monk,  and A. K. Parrot,  
``A dispersion analysis of finite element methods for Maxwell's equations," 
{ SIAM J. Sci. Stat. Comput.}, 
{ 15 (4)},   916 (1994).

\bibitem{dssy-ncell} J. Douglas, Jr.,  J. E. Santos,  D. Sheen and X. Ye,
`` Nonconforming  {Galerkin} methods based on quadrilateral elements for 
second order elliptic problems," 
{ RAIRO Mathematical Modelling and Numerical Analysis (M$2$AN) }
{33},  747 (1999).

\bibitem{santos_86} J. E. Santos,
``Elastic wave propagation in fluid-saturated porous media. Part I. The existence
and uniqueness results," {Mathematical Modelling and Numerical Analysis (M$2$AN) },
 { 20} (1),  113 (1986).
\bibitem{lions_69} J. L. Lions, 
`` Quelques methodes dde resolution des problems aux limites nonlineaires",
Dunod, Gauthier-Villars, Paris, (1969). 

\bibitem{santos_sinum_07} J. E. Santos and D. Sheen,
``Finite element methods for the simulation of waves in composite saturated
poroviscoelastic materials," { SIAM J. Numer. Anal.},
 { 45} (1),  389 (2007). 

\bibitem{zs07}F. I. Zyserman and J. E. Santos, 
``Analysis of the numerical dispersion of waves in  saturated poroelastic
 media," { Computer  Methods in Applied  Mechanics and Engineering}, 
{ 196},  4644 (2007).

\bibitem{dou00} J. Douglas, Jr., J. E. Santos, and D. Sheen,
``A nonconforming mixed finite element method for Maxwell's equations,"
{ Mathematical Models and Methods in Applied Sciences}, { 10},  593 (2000).

\bibitem{santos2dmt_98} J. E. Santos,
``Global and domain-decomposed mixed methods for the solution of
Maxwell's equations with application to magnetotellurics,"
 { Numer. Methods  Partial Differ. Eqs.}, 
{14},  1 (1998).

\bibitem{santos00} J. E. Santos and D. Sheen,
``On the existence and uniqueness of Maxwell's equations in bounded
domains with application to magnetotellurics,"
  { Mathematical Models and Methods in Applied Sciences}, 
 { 10},  615 (2000).

\bibitem{zys99} F. I. Zyserman, L. Guarracino and J. E. Santos,
``A hybridized mixed finite element domain decomposed
method for two-dimensional  magnetotelluric modelling," 
 { Earth, Planets and Space}, 
{ 51}, 297 (1999).

\bibitem{zys00} F. I. Zyserman and J. E. Santos,
``Parallel finite element algorithm for three--dimensional
magnetotelluric modelling," { J. Appl. Geophysics}, 
{ 44},   337 (2000).


\bibitem{johnston87} D. L. Johnston, J. Koplik and R. Dashen,
``Theory of dynamic permeability and tortuosity in fluid--saturated porous
media," {  J. Fluid Mechanics},  { 176}, 379 (1987).

\bibitem{gassmann51}
F. Gassmann, 
``Uber die elastizit\"at por\"oser medien," {  Vierteljahrsschrift 
der Naturforschenden Gessellschaft in Zurich}, { 96},  1 (1951).

\bibitem{santos92}
J. E. Santos, J.  Douglas  Jr, M. E.  Morley  and O. M.  Lovera,  
``Reflection and transmission coefficients in fluid-saturated porous media,"
 { J. Acoust. Soc. Amer.}, { 91},  1911 (1992).


\bibitem{biot56} M. A. Biot,
``Theory of deformation of a porous viscoelastic anisotropic solid,"
{ J. Appl. Physics}  27, 459 (1956).


\bibitem{biot62} M. A. Biot,  
``Mechanics of deformation and acoustic
propagation in porous media," { J. Appl. Physics} { 33 (4)},  1482 (1962).


\bibitem{liu76}  H. P. Liu and D. L. Anderson and H. Kanamori
``Velocity dispersion due to anelasticity; implications for seismology
and mantle composition," {  Geophys. J. R. Astr. Soc.} { l47},   
41 (1976).


\bibitem{seg_book_87} S. Ward and G. W. Hohmann, 
Electromagnetic theory for
geophysical applications, in 
Electromagnetic Methods in Applied Geophysics, Vol. 1, Theory, 
 M. N. Nabiguian, ed., { SEG Investigations in Geophysics} \# 3, (1987). 


\bibitem{santos_imajna} J. E. Santos,  J. Douglas, Jr., M. E. Morley, and O. M. Lovera,
``Finite element methods for a model for full waveform acoustic logging", 
 { IMA Journal of Numerical Analysis}, { 8}, 415 (1998).


\bibitem{gira} V. Girault  and P. Raviart,  Finite Element Methods for
Navier-Stokes Equations, {Springer--Verlag, Berlin}, 1986.


\bibitem{dsheenb} D. Sheen, 
``A generalized Green's theorem," { Appl. Math.
Lett.}, { 5},   95 (1992).

\bibitem{dsheenc} D. Sheen, 
``Approximation of electromagnetic fileds: Part I. Continuous problems," 
{SIAM J. Appl. Math.}, {6},  1716-1736 (1997).

\bibitem{buffa_01} A. Buffa and P. Ciarlet, 
``On traces for functional spaces related to maxwell's equations. II. Hodge
decompositions on the boundary of Lipschitiz polyhedra and applications," 
{Math. Methods Appl. Sci.}, {24 (1)},   31 (2001). 

\bibitem{buffa_02} A. Buffa and M. Costabel and D. Sheen, 
``On traces for  $H(\text{curl}, \O)$ in lipschitz domains," 
{J. math. Anal. Appl.}, {276 (2)},   845 (2002).

\bibitem{rav_05} C. L. Ravazzoli and J. E. Santos,   
``A theory for wave propagation in porous rocks saturated
by two-phase fluids  under variable pressure conditions," 
 { Bollettino di Geofisica teorica ed applicata}, 
{ 46 (4)}, 261 (2005).

\bibitem{ciarlet_78}, P. G. Ciarlet,
 The Finite Element Method for Elliptic Equations ,
{ North--Holland, Amsterdam }, 1978.

\bibitem{rduran}  R. Dur\'an, 
``Finite elements in $H(\text{curl})$ with
applications to elasticity," { Numer. Methods  Partial Differ.
Eqs.}, 4,   167 (1990).


\bibitem{dr} J. Douglas  Jr., and J. E. Roberts, 
``Global estimates for
mixed methods for second order elliptic equations," { Math. Comput.}, 44,  39 (1985).


\bibitem{tho} Thomas, J. M., 
``Sur l'analyse num\'erique des methodes
d'\'el\'ements finis hybrides et mixtes," Ph.D.\ thesis, Universit\'e
Pierre-et-Marie Curie, Paris, (1977).
 


\bibitem{duvaut-lions} G. Duvaut and J. L. Lions,
Les In{\'e}qualities en M{\'e}chanique et en Physique,
{ Dunod, Paris}, 1972.






\bibitem{nitsche81} J. A. Nitsche, 
``On Korn's second inequality," 
 { R. A. I. R. O Anal. Numer.},  { 15},  237 (1981).


\end{thebibliography}
\end{document}

