%\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}