%\documentclass[doublespacing]{elsart} \documentclass{elsart} \usepackage{subeqnarray,amsmath,amsfonts,graphicx} \usepackage{oldgerm} \usepackage{rotating} \bibliographystyle{elsart-num} \oddsidemargin -.5cm \journal{Comput. Methods Appl. Mech. Engrg.} \renewcommand\tilde{\widetilde} \renewcommand\a{\alpha} \def\b{\beta} \newcommand\bA{\mathbf{A}} \def\bmu{\boldsymbol\mu} \renewcommand\d{\delta} \newcommand\D{\Delta} \def\dfrac{\displaystyle\frac} \def\Sum{\displaystyle\sum} \def\ai{\left\langle\langle} \def\ad{\right\rangle\rangle} \def\b{\beta} \def\k{\kappa} \def\l{\ell} \def\D{\Delta} \def\d{\delta} \def\e{\varepsilon} \def\dint{\displaystyle\int} \def\g{\gamma} \def\G{\Gamma} \def\h{{\mathbf h}} \def\hf{ f} \def\hR{\widehat R} \def\hu{ u} \def\hz{\widehat z} \def\hZ{\widehat Z} \def\Im{\operatorname{Im}} \def\v{{\mathbf v}} \def\H{{\mathbf H}} \def\n{{\mathbf n}} \def\R{{\mathbb R}} \def\b{\beta} \def\lla{\langle\langle} \def\Lam{\Lambda} \def\lam{\lambda} \def\some{{x,y,z}} \def\Lip{\operatorname{Lip}} \def\tr{\operatorname{tr}} \def\curl{{\nabla\times}} \def\div{{\nabla\cdot}} \def\opcurl{\operatorname{curl}} \def\opdiv{\operatorname{div}} \def\nab{{\nabla}} \def\grad{{\nabla}} \def\g{\gamma} \def\q{\quad} \def\2q{\quad\quad} \def\3q{\quad\quad\quad} \def\s{\sigma} \def\sG{\mathcal{G}} \def\sL{\mathcal{L}} \def\sB{\mathcal{B}} \def\sA{\mathcal{A}} \def\cV{\mathcal{V}} \def\sm{\setminus} \def\NC{\mathcal{NC}} \def\sG{\mathcal{G}} \def\sm{\setminus} \def\t{\theta} \def\T{\Theta} \def\Tau{\mathcal T} \def\talpha{\widetilde\alpha} \def\tA{\tilde A} \def\tB{\tilde B} \def\tC{\tilde C} \def\V{\Vert} \def\z{\zeta} \def\xo{(x,\omega)} \def\<{\left\langle} \def\>{\right\rangle} \def\div{\nabla\cdot} \def\grad{\nabla} \def\lma{\left\langle} \def\la{\left\langle} \def\lla{\langle\langle} \def\rra{\rangle\rangle} \def\rma{\right\rangle} \def\ra{\right\rangle} \def\qq{\quad\quad} \def\dt{{\varDelta t}} \def\dx{{\varDelta x}} \def\dy{{\varDelta y}} \def\dz{{\varDelta z}} \def\BbbR{\mathbb R} \def\O{\Omega} \def\o{\omega} \def\Oh{{\mathcal O}} \def\mD{{\mathcal D}} \def\mS{{\mathcal S}} \def\mT{{\mathcal T}} \def\mC{{\mathcal C}} \def\mN{{\mathcal N}} \def\p{\partial} \def\for{\,\forall} \def\a{\alpha} \def\nb{\nonumber} \def\var{\varphi} \def\w{\omega} \def\G{\Gamma} \def\R{{\mathbb R}} \def\:{{\,:\,}} \def\Set#1{{\left\{#1\right\}}} \def\norm#1{{\left\|#1\right\|}} \def\qnorm#1{{\left|\left|\left|#1\right|\right|\right|}} \def\abs#1{{\left|#1\right|}} \def\th#1{\theta_{jk}^{#1}} \renewcommand{\hat}{\widehat} \newcommand{\pp}[2]{\frac{\partial {#1}}{\partial {\nu_{#2}}}} \newcommand{\pu}[2]{P_{\tau}u_{#1}^{#2}} \renewcommand{\Re}{\operatorname{Re}} \renewcommand{\Im}{\operatorname{Im}} % 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}{\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}\\ 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} u^f, \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*} \hskip4cm {\mathcal P} =\begin{pmatrix} \rho_b I &\rho_f I \\ \rho_f I& m I \end{pmatrix} \end{eqnarray*} where $I$ is the identity matrix in $R^2$, 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^s,v^s)_{\O_p} +(F^f,v^f)_{\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} $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*} \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 ${\bf M}$ is positive definite since it is associated with the strain energy density. Thus, 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^s = F^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}} = \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}\label{eq_20} &&C_1 = \text{min}\left(\sigma_{\text{min}} - \frac{L_0 \eta}{2 \kappa_{0,\text{max}}}, \frac{\eta}{\kappa_{0,\text{max}}}- \frac{L_0 \eta}{2 \kappa_{0,\text{max}}}\right) > 0. \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_3 \|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_3 \|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} Thus apply Gronwall's lemma in \eqref{eq_29} and note that all terms in the left-hand side of \eqref{eq_29} are nonnegative 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} Under the hypothesis made above on the positive definitess of the matrices ${\bf M}$, ${\mathcal D}$, and the validity of \eqref{eq_21}, 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} 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^s,\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 ${\mathcal A}_{\zeta,h}$ coercive and continuous on ${\mathcal Z}^h$. 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}$, ${\mathcal D}$, and the validity of \eqref{eq_21}, 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} Under the hypothesis of the positive definitess of the matrices ${\bf M}$, ${\mathcal D}$, and the validity of \eqref{eq_21}, 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))}\label{eq_102}\\ &&\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))} \nonumber\\ &&\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))} + \|(u^f-u^{f,h})\cdot \nu \|_{L^{2}(J,L^2(\G_p))} \\ &&\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))} +\|H_2\|^2_{L^{\infty}(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 u^f}{\p t}, \frac{\p u^f}{\p t} \right)_{\O_p} + \left\langle {\mathcal D} S_{\G_p}(\frac{\p u}{\p t}), S_{\G_p}(\frac{\p u^f}{\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{\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). \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} \end{proof} \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_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}