\documentclass[reqno]{amsart}
\usepackage{hyperref}
\usepackage{amssymb}
\usepackage{graphicx}

\AtBeginDocument{{\noindent\small
\emph{Electronic Journal of Differential Equations},
Vol. 2016 (2016), No. 184, pp. 1--22.\newline
ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu}
\thanks{\copyright 2016 Texas State University.}
\vspace{8mm}}

\begin{document}
\title[\hfilneg EJDE-2016/184\hfil Mathematical models of seismics]
{Mathematical models of seismics in composite media:
elastic and poro-elastic components}

\author[A. Meirmanov, M. Nurtas \hfil EJDE-2016/184\hfilneg]
{Anvarbek Meirmanov, Marat Nurtas}

\address{Anvarbek Meirmanov \hfill\break
School of Mathematical Sciences and Information Technology,
Yachay Tech,
Ibarra, Ecuador}
\email{ameirmanov@yachaytech.edu.ec}

\address{Marat Nurtas \hfill\break
School of mathematics and cibernetics,
Kazakh-British Technical University,
Almaty, Kazakhstan}
\email{marat\_nurtas@mail.ru}

\thanks{Submitted May 28, 2016. Published July 12, 2016.}
\subjclass[2010]{35B27, 46E35, 76R99}
\keywords{Acoustics; two-scale expansion method;
  full wave field inversion; 
\hfill\break\indent numerical simulation}

\begin{abstract}
 In the present paper we consider elastic and poroelastic media
 having a common interface. We derive the macroscopic mathematical
 models for seismic wave propagation through these two different
 media as a homogenization of the exact mathematical model at
 the microscopic level. They consist of seismic equations for
 each component and boundary conditions at the common interface,
 which separates different media. To do this we use the two-scale
 expansion method in the corresponding integral identities, defining
 the weak solution. We illustrate our results with the numerical
 implementations of the inverse problem for  the simplest model.
\end{abstract}

\maketitle
\numberwithin{equation}{section}
\allowdisplaybreaks


\section{Introduction}

This article is devoted to a  description of  seismic
wave propagation in composite media $Q\subset \mathbb{R}^{3}$,
consisting of the elastic medium $\Omega^{(0)}$, poroelastic medium
$\Omega$, which is perforated by a periodic system of pores filled
 with a fluid, and common interface $S^{(0)}$ between $\Omega^{(0)}$
 and $\Omega$ (see Figures  \ref{fig1}, \ref{fig2}). 
That is, $Q=\Omega\cup S^{(0)}\cup\Omega^{(0)}$
and $\Omega=\Omega_{f}\cup\Gamma\cup\Omega_s$, where $\Omega_s$
is a solid skeleton, $\Omega_{f}$ is a pore space (liquid domain),
and $\Gamma$ is a common boundary ``solid skeleton-liquid domain''.

The structure of the heterogeneous medium $Q$ is too complicated and
 makes hard a numerical simulation
of seismic waves propagation in multiscale media.  The main difficulty here
is a presence of both components (solid and liquid) in each sufficiently
small subdomain of $Q$. It requires to change
the governing equations (from Lame's equations to the Stokes equations)
at the scale of some tens microns.

There are two basic methods to describe physical processes in such media:
the phenomenological method and the asymptotical one which is based on
the upscaling approaches. The  phenomenological approach for waves
propagation through a poroelastic medium \cite{Berr1,Berr2}
leads, in particular, to Biot model \cite{Biot1}-\cite{Biot3}.
It based on the system of axioms (relations between the parameters of
the medium), which define the given physical process.  But, there can be
another system of axioms defining the same process.
Thus, it is necessary choose the correct  authenticity criterion of the
mathematical description of the process. It can be, for example, the
physical experiment. As a rule, each phenomenological model contains
some set of phenomenological constants. Therefore, one can achieve
agreement between the suggested theory and selected series of experiments
changing these parameters.

\begin{figure}[ht]
 \includegraphics[width=0.4\textwidth]{fig1} % 2222.jpg
\caption{Domain in consideration}\label{fig1}
\end{figure}

The second method, suggested by  Burridge and  Keller \cite{B-K} and
 Sanchez-Palencia \cite{S-P}, based on the homogenization.
It consists of:

(1) an exact description of the process at the microscopic level based
on the fundamental laws of continuum mechanics,and

(2) the rigorous homogenization of the obtained mathematical model.

To explain the method we consider a characteristic function
$\chi_0(\mathbf{x})$ of
the pore space $\Omega_{f}$. Let  $L$ is the characteristic size of the
physical domain in consideration,  $\tau$ is the characteristic time
of the physical process, $\rho^0$ is the
mean density of water,  and $g$ is acceleration due gravity.
In dimensionless variables
\[
\mathbf{x}\to\frac{\mathbf{x}}{L}, \quad
\mathbf{w}\to\alpha_{\tau}\frac{\mathbf{w}}{L}, \quad
t\to\frac{t}{\tau}, \quad
\mathbf{F}\to\frac{\mathbf{F}}{g},\quad
\rho\to\frac{\rho}{\rho^0},
\]
the dynamic system for the displacements $\mathbf{w}$ and pressure $p$
of the medium takes the form
\cite{B-K,S-P,AM1}:
\begin{gather}\label{1.1}
\varrho\frac{\partial^2\mathbf{w}}{\partial t^2}=
\nabla\cdot \mathbb{P}+\varrho\mathbf{F}, \\
\label{1.2}
\mathbb{P}=\chi_0\,\alpha_{\mu}\mathbb{D}
(x,\frac{\partial\mathbf{w}}{\partial t})+(1-\chi_0)\alpha_{\lambda}\,
\mathbb{D}(x,\mathbf{w})+\big(\chi_0\alpha_{\nu}(\nabla\,\cdot
\frac{\partial\mathbf{w}}{\partial t})-p\big) \mathbb{I}, \\
\label{1.3}
p+\alpha_{p}\nabla\cdot \mathbf{w}=0.
\end{gather}
Equations \eqref{1.1}-\eqref{1.3} are understood in the sense of distributions
as corresponding integral identities. They are equivalent to the Stokes equations
\begin{gather}\label{1.4}
\varrho_{f}\frac{\partial\mathbf{v}}{\partial t}=
\nabla\cdot \mathbb{P}_{f}+\varrho_{f}\mathbf{F},\quad
\frac{\partial p}{\partial t}+\alpha_{p,f}\nabla\cdot \mathbf{v}=0, \\
\label{1.5}
\mathbb{P}_{f}=\alpha_{\mu}\mathbb{D}(x,\mathbf{v})+
\big(\alpha_{\nu}(\nabla\,\cdot\mathbf{v})
-p\big)\mathbb{I}
\end{gather}
for the velocity $\mathbf{v}=\frac{\partial\mathbf{w}}{\partial t}$
and pressure $p$ in the pore space $\Omega_{f}$ and the Lame equations
\begin{gather}\label{1.6}
\varrho_s\frac{\partial^2\mathbf{w}}{\partial t^2}=
\nabla\cdot \mathbb{P}_s+\varrho_s\mathbf{F},\quad
p+\alpha_{p,s}\nabla\cdot \mathbf{w}=0, \\
\label{1.7}
\mathbb{P}_s=\alpha_{\lambda}\mathbb{D}(x,\mathbf{w})-p\mathbb{I}
\end{gather}
for the solid displacements $\mathbf{w}$ and pressure $p$  in $\Omega_s$.

At the common boundary $\Gamma$ velocities and normal tensions are continuous:
\begin{equation}\label{1.8}
\frac{\partial\mathbf{w}}{\partial t}=\mathbf{v},\quad
\mathbb{P}_s\cdot\mathbf{n}=\mathbb{P}_{f}\cdot\mathbf{n}.
\end{equation}
Here $\mathbf{n}$ is a unit normal to $\Gamma$.

In \eqref{1.1}-\eqref{1.8},
$\mathbb{D}(x,\mathbf{u})=\frac{1}{2}(\nabla\mathbf{u}+
\nabla\mathbf{u}^{*})$ is the symmetric part of
$\nabla\mathbf{u}$, $\mathbb{I}$ is a unit tensor,
$\mathbf{F}$ is a given vector of distributed mass forces,
\begin{gather*}
\alpha_{p}=\alpha_{p,f}\chi_0+\alpha_{p,s}(1-\chi_0),\quad
\varrho=\varrho_{f}\,\chi_0+\varrho_s(1-\chi_0),
\\
\alpha_{\tau}=\frac{ L}{g\tau^2}, \quad
\alpha_{\mu}=\frac{2\mu}{\alpha_{\tau}\tau Lg\,\rho^0},\quad
\alpha_{\lambda}=\frac{2\lambda}{\alpha_{\tau} Lg\,\rho^0},
\\
\alpha_{\nu}=\frac{2\nu}{\alpha_{\tau}\tau Lg\,\rho^0},\quad
\alpha_{p,f}=\frac{\varrho_{f} c_{f}^2}{\alpha_{\tau} Lg} ,\quad
\alpha_{p,s}=\frac{\varrho_s c_s^2}{\alpha_{\tau} Lg},
\end{gather*}
where $\mu$ is the dynamic viscosity, $\nu$ is the bulk viscosity,
$\lambda$ is  the elastic constant, $\varrho_{f}$ and $\varrho_s$
are the respective mean dimensionless densities of the liquid in
pores and the solid skeleton, correlated with the mean density of
water $\rho^0$, and $c_{f}$ and $c_s$ are the speed of
compression sound waves in the pore liquid and in the solid skeleton
respectively.

\begin{figure}[htbp]
\begin{center}
 \includegraphics[width=0.7\textwidth]{fig2} % 1.jpg}
\end{center}
\caption{The pore structure}\label{fig2}
\end{figure}

The  mathematical model \eqref{1.1}--\eqref{1.3} can not be
useful for practical needs, since the function $\chi_0$ changes
its value from 0 to 1 on the scale of a few microns. Fortunately,
the system possesses a natural small parameter
$\varepsilon=\frac{l}{L}$, where $l$  is the average
size of pores. Thus, the most suitable way to get a practically
significant mathematical model, which approximate
\eqref{1.1}--\eqref{1.3}, is a homogenization or upscaling.
That is, we suppose the $\varepsilon$-periodicity of the solid skeleton,
let $\varepsilon$ to be variable, and look for the limit in \eqref{1.1}--\eqref{1.3}
 as $\varepsilon\to 0$.

There are different homogenized (limiting)  systems, depending on of
 $\alpha_{\mu},\alpha_{\lambda}$, \dots
Some of these numbers might be small and some might be large. We may represent
 them as a power of
$\varepsilon$, or as functions depending on $\varepsilon$.

Let
\begin{gather*}
\mu_0=\lim_{\varepsilon\searrow 0} \alpha_{\mu}(\varepsilon), \quad
\nu_0=\lim_{\varepsilon\searrow 0} \alpha_{\nu}(\varepsilon), \quad
\lambda_0=\lim_{\varepsilon\searrow 0} \alpha_{\lambda}(\varepsilon), \\
c_{f,0}^2=\lim_{\varepsilon\searrow 0} \alpha_{p,f}(\varepsilon),\quad
c_{s,0}^2=\lim_{\varepsilon\searrow 0} \alpha_{p,s}(\varepsilon), \\
\mu_1=\lim_{\varepsilon\searrow 0}\frac{\alpha_{\mu}}{\varepsilon^2},\quad
\lambda_1=\lim_{\varepsilon\searrow 0}\frac{\alpha_{\lambda}}{\varepsilon^2}.
\end{gather*}
It is clear that the choice of these limits depend on our willing.
For example, for $\varepsilon=10^{-2}$ and
$\alpha=2\cdot10^{-1}$ we may state that $\alpha=2\cdot\varepsilon^{-\frac{1}{2}}$,
or $\alpha=0.02\cdot\varepsilon^0$. It is usual procedure when we neglect
some terms in differential equations with small coefficients and get more simple
equations, still describing the physical process.

The detailed analyses of all possible limiting regimes has been done
in \cite{AM1, AM2}.
To describe the seismic in two different media (elastic and
poroelastic),  having a common interface we must chose one of the two methods
discussed above. The first method suggests only some guesses, while
the second method has a clear algorithm for the derivation of the
boundary conditions. That is why we choose here the second method.

We derive new seismic
equations in each component (elastic and poroelastic) and the boundary
conditions on the common boundary. For these boundary conditions the very
little is known and only for the liquid filtration (see for example \cite{AnM}).

For three different sets of $\mu_0,\,\lambda_0,\dots $ for  each component
we derive three different mathematical models,
which describe the process with different degrees of approximation.

We start with the integral identities, defining the weak solution
$\mathbf{w}^{\varepsilon}$ and $p^{\varepsilon}$, and use the
two-scale expansion method \cite{B-L-P, B-P}, when we look for
the solution in the form
\begin{gather*}
\mathbf{w}^{\varepsilon}(\mathbf{x},t)=\mathbf{w}(\mathbf{x},t)+
\mathbf{W}_0(\mathbf{x},t,\frac{\mathbf{x}}{\varepsilon})+
\varepsilon\,\mathbf{W}_1(\mathbf{x},t,\frac{\mathbf{x}}{\varepsilon})
+o(\varepsilon), \\
p^{\varepsilon}(\mathbf{x},t)=p(\mathbf{x},t)+
P_0(\mathbf{x},t,\frac{\mathbf{x}}{\varepsilon})+
\varepsilon\,P_1(\mathbf{x},t,\frac{\mathbf{x}}{\varepsilon})+
o(\varepsilon)
\end{gather*}
with 1-periodic in the variable $\mathbf{y}$  functions
$\mathbf{W}_i(\mathbf{x},t,\mathbf{y})$,
$P_i(\mathbf{x},t,\mathbf{y})$, $i=0,1,\dots $

This method is rather heuristic and may lead to the wrong answer.
But our guesses are based upon the strong theory, suggested by G.
Nguetseng  \cite{GN1,GN2}. For the rigorous derivation of seismic
equations in poroelastic media, which dictate the correct two-scale
expansion, see \cite{AM1}.

Finally, to calculate limits as $\varepsilon \to 0$ in
corresponding integral identities, we apply the well-known result
\begin{equation}\label{1.9}
\lim_{\varepsilon\to 0}\int_{\Omega}F(\mathbf{x},
\frac{\mathbf{x}}{\varepsilon},t)\,dx\,dt=\int_{\Omega}\big(\int_{Y}
F(\mathbf{x},\mathbf{y},t)dy\big)\,dx\,dt
\end{equation}
for any smooth 1-periodic in the variable $\mathbf{y}\in Y$
function $F(\mathbf{x},\mathbf{y},t)$.

\section{Statement of the problem}

For the sake of simplicity we suppose that
$Q=\{\mathbf{x}=(x_1,x_2,x_3)\in \mathbb{R}^{3}:
x_3>0\}$, $\Omega^{(0)}=\{\mathbf{x}=(x_1,x_2,x_3)\in
\mathbb{R}^{3}: 0<x_3<H\}$,
$\Omega=\{\mathbf{x}=(x_1,x_2,x_3)\in \mathbb{R}^{3}:
x_3>H\}$, $\mathbf{F}=0$, and
\[
\alpha_{p,f}=\bar{c}_{f}^2,\,\,\alpha_{p,s}=\bar{c}_s^2.
\]

Let $Y$ be a unit cube in $\mathbb{R}^{3}$,
$Y=Y_{f}\cup\gamma\cup Y_s$. We assume that pore space
$\Omega^{\varepsilon}_{f}$ is a periodic repetition  in $\Omega$ of
the elementary cell $\varepsilon Y_{f}$, the solid skeleton
$\Omega^{\varepsilon}_s$ is a periodic repetition in $\Omega$ of
the elementary cell $\varepsilon Y_s$, and the boundary
$\Gamma^{\varepsilon}$ between a pore space and a solid skeleton is
a periodic repetition in $\Omega$ of the boundary $\varepsilon
\gamma$.
Detailed description of the sets  $Y_{f}$ and $Y_s$ is done in \cite{AM1}.
From these suppositions,
\[
\chi_0(\mathbf{x})=\chi^{\varepsilon}(\mathbf{x})=
\chi(\frac{\mathbf{x}}{\varepsilon}),
\]
where $\chi(\mathbf{y})$ is a 1-periodic function such that
$\chi(\mathbf{y})=1$ for $\mathbf{y}\in Y_{f}$ and $\chi(\mathbf{y})=0$
for $\mathbf{y}\in Y_s$.

For a fixed $\varepsilon>0$ the displacement vector $\mathbf{w}^{\varepsilon}$
and pressure $p^{\varepsilon}$ satisfy Lame's system
\begin{gather}\label{2.1}
\varrho^{(0)}_s\frac{\partial^2\mathbf{w}^{\varepsilon}}{\partial t^2}=
\nabla\cdot \mathbb{P}^{(0)}_s,\quad
p^{\varepsilon}+\bar{c}_{s,0}^2\nabla\cdot \mathbf{w}^{\varepsilon}=0,\\
\label{2.2}
\mathbb{P}^{(0)}_s=\alpha^{(0)}_{\lambda}\,
\mathbb{D}(x,\mathbf{w}^{\varepsilon})-p^{\varepsilon}\mathbb{I}
\end{gather}
in the domain $\Omega^{(0)}$ for $t>0$, and the system
\eqref{1.1}-\eqref{1.3}  with
$\chi_0=\chi^{\varepsilon}(\mathbf{x})$,
$\varrho=\varrho^{\varepsilon}=\varrho_{f}\,\chi^{\varepsilon}+
\varrho_s(1-\chi^{\varepsilon})$, and
$\alpha_{p}=\alpha^{\varepsilon}_{p}=\alpha_{p,f}\chi^{\varepsilon}+
\alpha_{p,s}(1-\chi^{\varepsilon})$ in the domain $\Omega$ for
$t>0$.

On the common boundary
$S^{(0)}=\{\mathbf{x}=(x_1,x_2,x_3)\in \mathbb{R}^{3}:
x_3=H\}$  the displacement vector and normal tensions are
continuous:
\begin{gather}\label{2.3}
\lim_{\substack{\mathbf{x}\to \mathbf{x}^0\\
\mathbf{x}\in \Omega^{(0)}}}
\mathbf{w}^{\varepsilon}(\mathbf{x},t)=
\lim_{\substack{\mathbf{x}\to \mathbf{x}^0\\
\mathbf{x}\in \Omega}}\mathbf{w}^{\varepsilon}(\mathbf{x},t),
\quad \mathbf{x}^0\in S^{(0)}, \\
\label{2.4}
\lim_{\substack{\mathbf{x}\to \mathbf{x}^0\\
\mathbf{x}\in \Omega^{(0)}}}\mathbb{P}^{(0)}(\mathbf{x},t)
\cdot\mathbf{e}_3=
\lim_{\substack{\mathbf{x}\to \mathbf{x}^0\\
\mathbf{x}\in \Omega}}\mathbb{P}(\mathbf{x},t)
\cdot\mathbf{e}_3,\quad,\mathbf{x}^0\in S^{(0)},
\end{gather}
where $\mathbf{e}_3=(0,0,1)$.

The problem is complemented with the boundary condition
\begin{equation}\label{2.5}
\mathbb{P}^{(0)}_s\cdot\mathbf{e}_3=
-p^0(\mathbf{x}',t)\mathbf{e}_3,\quad \mathbf{x}'=(x_1,x_2)
\end{equation}
on the outer boundary $S=\{\mathbf{x}=(x_1,x_2,x_3)\in
\mathbb{R}^{3}: x_3=0\}$  for $t>0$ and homogeneous initial
conditions
\begin{equation}\label{2.6}
\mathbf{w}^{\varepsilon}(\mathbf{x},0)=
\frac{\partial\mathbf{w}^{\varepsilon}}{\partial t}(\mathbf{x},0)=0.
\end{equation}
Let $\varsigma(\mathbf{x})$ be the characteristic function of the domain $\Omega$ and
\[
\tilde{\varrho}^{\varepsilon}=(1-\varsigma)\varrho^{(0)}_s+
\varsigma\varrho^{\varepsilon},\,\,\tilde{\alpha}_{p}^{\varepsilon}=
(1-\varsigma)\bar{c}_{s,0}^2+\varsigma\,\alpha_{p}^{\varepsilon}.
\]
Then the above formulated problem takes the form
\begin{gather}\label{2.7}
\tilde{\varrho}^{\varepsilon}
\frac{\partial^2\mathbf{w}^{\varepsilon}}{\partial t^2}=
\nabla\cdot \big((1-\varsigma)\mathbb{P}^{(0)}_s+\varsigma\mathbb{P}\big), \\
\label{2.8}
p^{\varepsilon}+\tilde{\alpha}_{p}^{\varepsilon}
\nabla\cdot \mathbf{w}^{\varepsilon}=0, \\
\label{2.9}
\mathbb{P}=\chi^{\varepsilon}\,\alpha_{\mu}\mathbb{D}
(x,\frac{\partial\mathbf{w}^{\varepsilon}}{\partial t})+
(1-\chi^{\varepsilon})\alpha_{\lambda}\,
\mathbb{D}(x,\mathbf{w}^{\varepsilon})-
\big(\chi^{\varepsilon}\frac{\alpha_{\nu}}{\bar{c}_{f}^2}
\frac{\partial p^{\varepsilon}}{\partial t}+p^{\varepsilon}\big)\mathbb{I},
\end{gather}
where in \eqref{2.9} we have used the consequence of \eqref{2.8} in the form
\[
\varsigma\,\chi^{\varepsilon}\alpha_{\nu}(\nabla\,\cdot
\frac{\partial\mathbf{w}^{\varepsilon}}{\partial t})=
-\varsigma\,\chi^{\varepsilon}\frac{\alpha_{\nu}}{\bar{c}_{f}^2}
\frac{\partial p^{\varepsilon}}{\partial t}.
\]
Equation \eqref{2.7} is understood in the sense of distributions.
That is, for any smooth functions $\mathbf{\varphi}$ with a
compact support in $\overline{Q}$ the following integral identity
\begin{equation}\label{2.10}
\int_{Q_{T}}\Big(\tilde{\varrho}^{\varepsilon}
\frac{\partial^2\mathbf{w}^{\varepsilon}}{\partial t^2}\cdot\mathbf{\varphi}+
\big((1-\varsigma)\mathbb{P}^{(0)}_s+\varsigma\mathbb{P}\big):
\mathbb{D}(x,\mathbf{\varphi})+\nabla\cdot(p^0\mathbf{\varphi})\Big)\,dx\,dt=0
\end{equation}
holds. We call such solution the \emph{weak solution}.

In \eqref{2.10} $Q_{T}=Q\times(0,T)$ and the convolution
$\mathbb{A}:\mathbb{B}$ of  two tensors $\mathbb{A}=(A_{ij})$ and
$\mathbb{B}=(B_{ij})$ is defined as
$\mathbb{A}:\mathbb{B}=\mbox{tr}(\mathbb{A}\cdot\mathbb{B})=
\sum_{i,j=1}^{3}A_{ij}B_{ji}$.

Using standard methods one can prove that for any positive
$\varepsilon >0$  and given smooth function $p^0$ there exists a
unique weak solution of the problem \eqref{2.7}-\eqref{2.9} which
makes sense to the integral identity \eqref{2.10}.
We look for the limit of the weak solutions  for the following cases:
\begin{itemize}
\item[(I)] $\mu_0=\lambda_0=\lambda^{(0)}_0=0$,
$\mu_1=\lambda_1=\infty$, $0\leqslant\nu_0<\infty$,\
$\lambda^{(0)}_0=\lim_{\varepsilon\to0}\alpha^{(0)}_{\lambda}$;

\item[(II)] $\mu_0=\lambda_0=\lambda^{(0)}_0=\mu_1=\nu_0=0$,
$\lambda_1=\infty$;

\item[(III)] $\mu_0=\nu_0=0$, $0<\lambda_0$, $\lambda^{(0)}_0$, $\mu_1<\infty$.
\end{itemize}

\section{Homogenization: case (I)}

According to \cite{AM2}, the two-scale expansion for the weak
solution of  the problem \eqref{2.7}-\eqref{2.9} under conditions
$(\textbf{I})$ has the form
\begin{equation}\label{3.1}
\mathbf{w}^{\varepsilon}(\mathbf{x},t)=
\mathbf{w}(\mathbf{x},t)+o(\varepsilon),\quad
p^{\varepsilon}(\mathbf{x},t)=p(\mathbf{x},t)+o(\varepsilon),
\end{equation}
where $\lim_{\varepsilon\to 0}o(\varepsilon)=0$.

The substitution \eqref{3.1} into \eqref{2.10} results in the integral identity
\begin{equation} \label{3.2}
\begin{aligned}
&\int_{\Omega_{T}}\Big(\big(\chi(\frac{\mathbf{x}}{\varepsilon})\varrho_{f}+
\big(1-\chi(\frac{\mathbf{x}}{\varepsilon})\big)\varrho_s\big)
\frac{\partial^2\mathbf{w}}{\partial t^2}\cdot\mathbf{\varphi}-
(\chi(\frac{\mathbf{x}}{\varepsilon})\frac{\alpha_{\nu}}{\bar{c}_{f}^2}
\frac{\partial p}{\partial t}+p)\nabla\cdot \mathbf{\varphi}\Big)\,dx\,dt\\
&+\int_{Q_{T}}\nabla\cdot(p^0\mathbf{\varphi})\,dx\,dt+
\int_{\Omega^{(0)}_{T}}\big(\varrho^{(0)}_s
\frac{\partial^2\mathbf{w}}{\partial t^2}\cdot\mathbf{\varphi}-
p(\nabla\cdot \mathbf{\varphi})\big)\,dx\,dt=o(\varepsilon).
\end{aligned}
\end{equation}
Now we use \eqref{1.9} and after the limit in \eqref{3.2} as
$\varepsilon\to 0$   arrive at the  integral identity
\begin{equation} \label{3.3}
\begin{aligned}
&\int_{\Omega_{T}}\big(\hat{\varrho}\frac{\partial^2\mathbf{w}}{\partial t^2}
\cdot\mathbf{\varphi}-(m\frac{\nu_0}{\bar{c}_{f}^2}
\frac{\partial p}{\partial t}+p)\nabla\cdot \mathbf{\varphi}\big)\,dx\,dt+
\int_{Q_{T}}\nabla\cdot(p^0\mathbf{\varphi})\,dx\,dt\\
&+ \int_{\Omega^{(0)}_{T}}\big(\varrho^{(0)}_s
\frac{\partial^2\mathbf{w}}{\partial t^2}\cdot\mathbf{\varphi}-
p(\nabla\cdot \mathbf{\varphi})\big)\,dx\,dt=0,
\end{aligned}
\end{equation}
where $\hat{\varrho}=m\varrho_{f}+(1-m)\varrho_s$ and
 $ m=\int_{Y}\chi(\mathbf{y})dy$.

Next we rewrite \eqref{2.8} as
\begin{equation}\label{3.4}
\big(\frac{(1-\varsigma)}{\bar{c}_{s,0}^2}+
\frac{\varsigma\,\chi(\frac{\mathbf{x}}{\varepsilon})}{\bar{c}_{f}^2}+
\frac{\varsigma\big(1-\chi(\frac{\mathbf{x}}{\varepsilon})\big)}
{\bar{c}_s^2}\big) p^{\varepsilon}+
\nabla\cdot \mathbf{w}^{\varepsilon}=0\,.
\end{equation}
We multiply the result by a smooth function $\psi(\mathbf{x},t)$
with a compact support in $Q$, and integrate by parts over domain
$Q_{T}$:
\begin{equation}\label{3.5}
\int_{Q_{T}}\Big(\psi\big(\frac{(1-\varsigma)}{\bar{c}_{s,0}^2}+
\frac{\varsigma\,\chi(\frac{\mathbf{x}}{\varepsilon})}{\bar{c}_{f}^2}+
\frac{\varsigma\big(1-\chi(\frac{\mathbf{x}}{\varepsilon})\big)}
{\bar{c}_s^2}\big) p^{\varepsilon}-
\nabla\,\psi\cdot \mathbf{w}^{\varepsilon}\Big)\,dx\,dt=0.
\end{equation}
As above, we substitute \eqref{3.1} into \eqref{3.5} and pass to the limit as
$\varepsilon\to 0$:
\begin{equation}\label{3.6}
\int_{Q_{T}}\Big(\psi\big(\frac{(1-\varsigma)}{\bar{c}_{s,0}^2}+
\frac{\varsigma\,m}{\bar{c}_{f}^2}+
\frac{\varsigma(1-m)}{\bar{c}_s^2}\big)
p-\nabla\,\psi\cdot \mathbf{w}\Big)\,dx\,dt=0.
\end{equation}
Integral identities \eqref{3.3} and \eqref{3.6},  complemented with initial
 conditions
\begin{equation}\label{3.7}
\mathbf{w}(\mathbf{x},0)=
\frac{\partial\mathbf{w}}{\partial t}(\mathbf{x},0)=0,
\end{equation}
form mathematical model (I) of seismic in composite media.

In fact, these identities contain the differential equations in
$\Omega$  and $\Omega^{(0)}$ and the boundary conditions on $S$ and
$S^{(0)}$.

Let $\mathbf{\varphi}$ be a smooth function with a compact
support  in $\Omega^{(0)}$. Rewriting \eqref{3.3} as
\[
\int_{\Omega^{(0)}_{T}}(\varrho^{(0)}_s
\frac{\partial^2\mathbf{w}}{\partial t^2}+
\nabla p)\cdot\mathbf{\varphi}\,dx\,dt=0
\]
and using the arbitrary choice of $\mathbf{\varphi}$ we conclude that
\begin{equation}\label{3.8}
\varrho^{(0)}_s\frac{\partial^2\mathbf{w}}{\partial t^2}=-\nabla p
\end{equation}
in the domain $\Omega^{(0)}$ for $t>0$.

For functions $\mathbf{\varphi}$ with a compact support in $\Omega$,
 \eqref{3.3} implies
\begin{equation}\label{3.9}
\hat{\varrho}\frac{\partial^2\mathbf{w}}{\partial t^2}=
-\nabla(p+m\frac{\nu_0}{\bar{c}_{f}^2}\frac{\partial p}{\partial t}),\quad
\hat{\varrho}=m\varrho_{f}+(1-m)\varrho_s
\end{equation}
in the domain $\Omega$ for $t>0$.

Now, if we choose
$\mathbf{\varphi}=(\varphi_1,\varphi_2,\varphi_3)$ with  a
compact support in $Q$ and $\varphi_3(\mathbf{x},t)\neq 0$ for
$\mathbf{x}\in S^{(0)}$, then the integration by parts in
\eqref{3.3} together with \eqref{3.8} and \eqref{3.9} result in
\[
\int_{S^{(0)}_{T}}\big(p^{-}-(p^{+}+m\frac{\nu_0}{\bar{c}_{f}^2}
\frac{\partial p^{+}}{\partial t})\big)\varphi_3\,dx\,dt=0,
\]
where
\[
p^{-}(x_1,x_2,t)=p(x_1,x_2,H-0,t),\quad
p^{+}(x_1,x_2,t)=p(x_1,x_2,H+0,t).
\]
Therefore,
\begin{equation}\label{3.10}
\lim_{\substack{\mathbf{x}\to \mathbf{x}^0\\
\mathbf{x}\in \Omega^{(0)}}}
p(\mathbf{x},t)=\lim_{\substack{
\mathbf{x}\to \mathbf{x}^0\\
\mathbf{x}\in \Omega}}\big(p(\mathbf{x},t)+
m\frac{\nu_0}{\bar{c}_{f}^2}
\frac{\partial p}{\partial t}(\mathbf{x},t)\big),
\,\,\,\mathbf{x}^0\in S^{(0)}.
\end{equation}
Finally, for functions $\mathbf{\varphi}$ with a compact support
in $\Omega^{(0)}$ and $\varphi_3(\mathbf{x},t)\neq 0$ for
$\mathbf{x}\in S$ the integration by parts in \eqref{3.3}
together with \eqref{3.8} result in
\[
\int_{S_{T}}(p-p^0)\varphi_3\,dx\,dt=0,
\]
or
\begin{equation}\label{3.11}
p(\mathbf{x},t)=p^0(\mathbf{x},t),\,\,\,\mathbf{x}\in S.
\end{equation}
In the same way as above, it can be shown that \eqref{3.6}
implies continuity equations
\begin{equation}\label{3.12}
\frac{1}{\bar{c}_{s,0}^2} p+\nabla\cdot \mathbf{w}=0
\end{equation}
and
\begin{equation}\label{3.13}
\big(\frac{m}{\bar{c}_{f}^2}+\frac{(1-m)}{\bar{c}_s^2}\big)\,
p+\nabla\cdot \mathbf{w}=0
\end{equation}
in the domains $\Omega^{(0)}$ and $\Omega$ respectively, and the boundary
condition
\begin{equation}\label{3.14}
\lim_{\substack{\mathbf{x}\to \mathbf{x}^0\\
\mathbf{x}\in \Omega^{(0)}}}
\mathbf{e}_3\cdot\mathbf{w}(\mathbf{x},t)=
\lim_{\substack{\mathbf{x}\to \mathbf{x}^0\\
\mathbf{x}\in \Omega}}
\mathbf{e}_3\cdot\mathbf{w}(\mathbf{x},t),
\quad \mathbf{x}^0\in S^{(0)}
\end{equation}
on the common boundary $S^{(0)}$.

Differential equations \eqref{3.8}, \eqref{3.9}, \eqref{3.12},  and
\eqref{3.13}, boundary conditions \eqref{3.10}, \eqref{3.11}, and
\eqref{3.14}, and initial conditions \eqref{3.7} constitute the
mathematical model (I) in its differential form.

\section{Homogenization: case (II)}

For this case we put
\begin{equation} \label{4.1}
\begin{gathered}
\mathbf{w}^{\varepsilon}(\mathbf{x},t)=
(1-\varsigma)\mathbf{w}(\mathbf{x},t)+
\varsigma\chi(\frac{\mathbf{x}}{\varepsilon})\mathbf{w}^{(f,\varepsilon)}
(\mathbf{x},t)+
\varsigma\big(1-\chi(\frac{\mathbf{x}}{\varepsilon})\big)
\mathbf{w}_s(\mathbf{x},t)+o(\varepsilon),\\
p^{\varepsilon}(\mathbf{x},t)=p(\mathbf{x},t)+o(\varepsilon),
\end{gathered}
\end{equation}
where
\[
\mathbf{w}^{(f,\varepsilon)}(\mathbf{x},t)=
\mathbf{W}^{(f)}(\mathbf{x},t,\frac{\mathbf{x}}{\varepsilon}),
\]
and $\mathbf{W}^{(f)}(\mathbf{x},t,\mathbf{y})$ is a
1-periodic  in the variable $\mathbf{y}$ function.

The substitution \eqref{4.1} into \eqref{2.10} results in the integral identity
\begin{equation} \label{4.2}
\begin{aligned}
&\int_{Q_{T}}\nabla\cdot(p^0\mathbf{\varphi})\,dx\,dt+
\int_{\Omega^{(0)}_{T}}\big(\varrho^{(0)}_s
\frac{\partial^2\mathbf{w}}{\partial t^2}\cdot\mathbf{\varphi}-
p(\nabla\cdot \mathbf{\varphi})\big)\,dx\,dt\\
&+\int_{\Omega_{T}}\Big(\big(\varrho_{f}\chi(\frac{\mathbf{x}}{\varepsilon})
\frac{\partial^2\mathbf{W}^{(f)}}{\partial t^2}
(\mathbf{x},t,\frac{\mathbf{x}}{\varepsilon})+
\varrho_s\big(1-\chi(\frac{\mathbf{x}}{\varepsilon})\big)
\frac{\partial^2\mathbf{w}_s}{\partial t^2}\big)\cdot\mathbf{\varphi}-
p(\nabla\cdot \mathbf{\varphi})\Big)\,dx\,dt\\
&= -\int_{\Omega_{T}}\alpha_{\mu}\,\chi(\frac{\mathbf{x}}{\varepsilon})\mathbb{D}
(x,\frac{\partial\mathbf{w}^{(f,\varepsilon)}}{\partial t}):
\mathbb{D}(x,\mathbf{\varphi})\,dx\,dt+o(\varepsilon),
\end{aligned}
\end{equation}
which holds  for any smooth function $\mathbf{\varphi}(\mathbf{x},t)$.
Let
\[
\mathbf{w}^{(f)}(\mathbf{x},t)
=\varsigma\,\lim_{\varepsilon\to 0}\chi(\frac{\mathbf{x}}{\varepsilon})
\mathbf{W}^{(f)}(\mathbf{x},t,\frac{\mathbf{x}}{\varepsilon})=
\varsigma\,\int_{Y_{f}}\mathbf{W}^{(f)}(\mathbf{x},t,\mathbf{y})dy
\]
be the weak limit of the sequence
$\{\mathbf{w}^{\varepsilon}\}$.  Then after the limit as
$\varepsilon \to 0$ we arrive at the integral identity
\begin{equation} \label{4.3}
\begin{aligned}
&\int_{Q_{T}}\nabla\cdot(p^0\mathbf{\varphi})\,dx\,dt+
\int_{\Omega^{(0)}_{T}}\big(\varrho^{(0)}_s
\frac{\partial^2\mathbf{w}}{\partial t^2}\cdot\mathbf{\varphi}-
p(\nabla\cdot \mathbf{\varphi})\big)\,dx\,dt\\
&=\int_{\Omega_{T}}\Big(\big(\varrho_{f}
\frac{\partial^2\mathbf{w}^{(f)}}{\partial t^2}
+\varrho_s(1-m)\frac{\partial^2\mathbf{w}_s}{\partial t^2}\big)
\cdot\mathbf{\varphi}-p(\nabla\cdot \mathbf{\varphi})\Big)\,dx\,dt=0.
\end{aligned}
\end{equation}
Note that the term $\alpha_{\mu}\mathbb{D}
(x,\frac{\partial\mathbf{w}^{(f,\varepsilon)}}{\partial t})$  in
the right-hand side of \eqref{4.2} converges to zero due to the
supposition $\lim_{\varepsilon\searrow 0} \alpha_{\mu}=
\lim_{\varepsilon\searrow 0}\frac{\alpha_{\mu}}{\varepsilon}=0$:
\[
\alpha_{\mu}\mathbb{D}\big(x,\frac{\partial\mathbf{w}^{(f,\varepsilon)}}{\partial t}
(\mathbf{x},t)\big)=\alpha_{\mu}\mathbb{D}
\big(x,\frac{\partial\mathbf{W}^{(f)}}{\partial t}
(\mathbf{x},t,\frac{\mathbf{x}}{\varepsilon})\big)+
\frac{\alpha_{\mu}}{\varepsilon}\mathbb{D}
\big(y,\frac{\partial\mathbf{W}^{(f)}}{\partial t}
(\mathbf{x},t,\frac{\mathbf{x}}{\varepsilon})\big).
\]
The substitution of \eqref{4.1} into the continuity equation
\eqref{2.8}  leads to the integral identity
\begin{equation} \label{4.4}
\begin{aligned}
&\int_{Q_{T}}\Big(\psi\big(\frac{(1-\varsigma)}{\bar{c}_{s,0}^2}+
\frac{\varsigma\,\chi^{\varepsilon}}{\bar{c}_{f}^2}+
\frac{\varsigma(1-\chi^{\varepsilon})}{\bar{c}_s^2}\big) p\\
&-\nabla \psi\cdot\big((1-\varsigma)\mathbf{w}+
\varsigma \chi^{\varepsilon}\mathbf{W}^{(f)}
(\mathbf{x},t,\frac{\mathbf{x}}{\varepsilon})+
\varsigma(1-\chi^{\varepsilon})\mathbf{w}_s\big)\Big)\,dx\,dt=o(\varepsilon).
\end{aligned}
\end{equation}
The limit here as $\varepsilon \to 0$ results in
\begin{equation} \label{4.5}
\begin{aligned}
\int_{Q_{T}}\Big(\psi\big(\frac{(1-\varsigma)}{\bar{c}_{s,0}^2}+
\frac{\varsigma\,m}{\bar{c}_{f}^2}+
\frac{\varsigma(1-m)}{\bar{c}_s^2}\big) p+\\
\nabla\,\psi\cdot\big((1-\varsigma)\mathbf{w}+
\varsigma\,\mathbf{w}^{(f)}+\varsigma(1-m)\mathbf{w}_s\big)\Big)\,dx\,dt=0.
\end{aligned}
\end{equation}
As in previous section we conclude that integral identities
\eqref{4.3} and \eqref{4.5}  imply differential equations in
$\Omega^{(0)}$ and $\Omega$ and boundary conditions on the
boundaries $S^{(0)}$ and $S$.

Namely, in the domain $\Omega^{(0)}$ the displacements vector
$\mathbf{w}$ and  pressure $p$ of the solid component satisfy
the seismic system
\begin{gather}\label{4.6}
\varrho^{(0)}_s\frac{\partial^2\mathbf{w}}{\partial t^2}=-\nabla p, \\
\label{4.7}
\frac{1}{\bar{c}_{s,0}^2} p+\nabla\cdot \mathbf{w}=0.
\end{gather}
In the domain $\Omega$ the displacements vector
$\mathbf{w}_s$ of the solid  component, displacements vector
$\mathbf{w}^{(f)}$ of the liquid component, and pressure $p$ of
the medium satisfy the seismic system
\begin{gather}\label{4.8}
\varrho_{f}\frac{\partial^2\mathbf{w}^{(f)}}{\partial t^2}
+\varrho_s(1-m)\frac{\partial^2\mathbf{w}_s}{\partial t^2}=-\nabla p, \\
\label{4.9}
\big(\frac{m}{\bar{c}_{f}^2}+\frac{(1-m)}{\bar{c}_s^2}\big) p+
\nabla\cdot\big(\mathbf{w}^{(f)}+(1-m)\mathbf{w}_s\big)=0.
\end{gather}
On the common boundary $S^{(0)}$ the displacements vectors
$\mathbf{w}$,  $\mathbf{w}_s$, and $\mathbf{w}^{(f)}$
and pressure $p$ satisfy continuity conditions
\begin{gather} \label{4.10}
\lim_{\substack{\mathbf{x}\to \mathbf{x}^0\\
\mathbf{x}\in \Omega^{(0)}}}
p(\mathbf{x},t)=\lim_{\substack{
\mathbf{x}\to \mathbf{x}^0\\
\mathbf{x}\in \Omega}}p(\mathbf{x},t),
\\
\label{4.11}
\lim_{\substack{\mathbf{x}\to \mathbf{x}^0\\
\mathbf{x}\in \Omega^{(0)}}}
\mathbf{e}_3\cdot\mathbf{w}(\mathbf{x},t)=
\lim_{\substack{\mathbf{x}\to \mathbf{x}^0\\
\mathbf{x}\in \Omega}}
\mathbf{e}_3\cdot\big(\mathbf{w}^{(f)}(\mathbf{x},t)+
(1-m)\mathbf{w}_s(\mathbf{x},t)\big).
\end{gather}
Finally, on the outer boundary $S$,
\begin{equation}\label{4.12}
p(\mathbf{x},t)=p^0(\mathbf{x},t).
\end{equation}
As above, we have to add the initial conditions:
\begin{equation} \label{4.13}
\begin{aligned}
\mathbf{w}(\mathbf{x},0)
&=\frac{\partial\mathbf{w}}{\partial t}(\mathbf{x},0)
 =\mathbf{w}_s(\mathbf{x},0)
 =\frac{\partial\mathbf{w}_s}{\partial t}(\mathbf{x},0)\\
&=\mathbf{w}^{(f)}(\mathbf{x},0)=
\frac{\partial\mathbf{w}^{(f)}}{\partial t}(\mathbf{x},0)=0.
\end{aligned}
\end{equation}
The obtained system of differential equations and boundary and
initial conditions  is still incomplete. We have no differential
equation for the liquid displacements $\mathbf{w}^{(f)}$. To
find the missing equation we pass to the limit
$\varepsilon\to 0$ in \eqref{4.2} with test functions
$\mathbf{\varphi}^{\varepsilon}$ in the form
\[
\mathbf{\varphi}^{\varepsilon}(\mathbf{x},t)=h(\mathbf{x},t)
\mathbf{\varphi}_0(\frac{\mathbf{x}}{\varepsilon}),
\]
where $h$ is a smooth function with a compact support in $\Omega$
and  $\mathbf{\varphi}_0(\mathbf{y})$ is a smooth function
with a compact support in $Y_{f}$ (that is
$\mathbf{\varphi}^{\varepsilon}$ vanishes outside of the pore
space $\Omega_{f}$).

For an arbitrary function $\mathbf{\varphi}_0(\mathbf{y})$
the term $p\nabla\cdot \mathbf{\varphi}^{\varepsilon}$
becomes unbounded as $\varepsilon\to 0$:
\[
\nabla\cdot \mathbf{\varphi}^{\varepsilon}=\big(\nabla_{x}\,h(\mathbf{x},t)\big)
\cdot \mathbf{\varphi}_0(\frac{\mathbf{x}}{\varepsilon})+
\frac{1}{\varepsilon}h(\mathbf{x},t)\big(\nabla_{y}\cdot \mathbf{\varphi}
_0(\frac{\mathbf{x}}{\varepsilon})\big).
\]
Therefore, we require that conditions
\begin{gather}\label{4.14}
\nabla_{y}\cdot \mathbf{\varphi}_0(\mathbf{y})=0,\quad \mathbf{y}\in Y_{f}, \\
\label{4.15}
\mathbf{\varphi}_0(\mathbf{y})=0,\quad \mathbf{y}\in \gamma
\end{gather}
hold.

The term $\alpha_{\mu}\mathbb{D}
(x,\frac{\partial\mathbf{w}^{(f,\varepsilon)}}{\partial t}):
\mathbb{D}(x,\mathbf{\varphi}^{\varepsilon})$ in the right-hand
side of \eqref{4.2}  converges to zero because of the  assumptions
$\lim_{\varepsilon\searrow 0}
\alpha_{\mu}=\lim_{\varepsilon\searrow
0}\frac{\alpha_{\mu}}{\varepsilon} =\lim_{\varepsilon\searrow
0}\frac{\alpha_{\mu}}{\varepsilon^2}=0$:
\begin{align*}
&\alpha_{\mu}\mathbb{D}\big(x,\frac{\partial\mathbf{w}
^{(f,\varepsilon)}}{\partial t}
(\mathbf{x},t)\big):\mathbb{D}(x,\mathbf{\varphi}^{\varepsilon})\\
&=\alpha_{\mu}\Big(\mathbb{D}\big(x,\frac{\partial\mathbf{W}^{(f)}}{\partial t}
(\mathbf{x},t,\frac{\mathbf{x}}{\varepsilon})\big)+
\frac{1}{\varepsilon}\mathbb{D}\big(y,\frac{\partial\mathbf{W}^{(f)}}{\partial t}
(\mathbf{x},t,\frac{\mathbf{x}}{\varepsilon})\big)\Big):\\
&\quad \Big(\frac{1}{2}\big((\nabla\,h)\otimes\mathbf{\varphi}_0+
\mathbf{\varphi}_0\otimes(\nabla\,h)\big)+\frac{h}{\varepsilon}\mathbb{D}
\big(y,\mathbf{\varphi}_0(\frac{\mathbf{x}}{\varepsilon})\big)\Big)=
o(\varepsilon).
\end{align*}
Here a matrix $\mathbf{a}\otimes\mathbf{b}$ is defined as
\[
(\mathbf{a}\otimes\mathbf{b})\cdot\mathbf{c}=\mathbf{a}
(\mathbf{b}\cdot\mathbf{c}),
\]
for any vectors $\mathbf{a}$, $\mathbf{b}$, and $\mathbf{c}$.

Thus, the limit as $\varepsilon\to 0$ in \eqref{4.2} results in the integral identity
\begin{equation} \label{4.16}
\begin{aligned}
&\int_{\Omega_{T}}\Big(\int_{Y_{f}}\big(\varrho_{f}
\frac{\partial^2\mathbf{W}^{(f)}}{\partial t^2}\,h-p\,
(\nabla\cdot h)\big)\cdot\mathbf{\varphi}_0dy\Big)\,dx\,dt\\
&=\int_{\Omega_{T}}h(\mathbf{x},t)\Big(\int_{Y_{f}}\big(\varrho_{f}
\frac{\partial^2\mathbf{W}^{(f)}}{\partial t^2}+\nabla p\big)\cdot
\mathbf{\varphi}_0dy\Big)\,dx\,dt=0,
\end{aligned}
\end{equation}
which holds for any smooth function $h(\mathbf{x},t)$ with
a compact  support in $\Omega$ and for any smooth solenoidal
function $\mathbf{\varphi}_0(\mathbf{y})$ with a compact
support in $Y_{f}$.

By arbitrary choice of $h(\mathbf{x},t)$, \eqref{4.16} implies
\begin{equation}\label{4.17}
\int_{Y_{f}}\big(\varrho_{f}
\frac{\partial^2\mathbf{W}^{(f)}}{\partial t^2}+
\nabla p\big)\cdot\mathbf{\varphi}_0dy=0.
\end{equation}
This identity means that the function
$\big(\varrho_{f}
\frac{\partial^2\mathbf{W}^{(f)}}{\partial
t^2}+\nabla p\big)$ is orthogonal  to any solenoidal function.
Therefore there exists some 1-periodic in the variable
$\mathbf{y}$ function $\Pi(\mathbf{x},t,\mathbf{y})$
such that
\begin{equation}\label{4.18}
\varrho_{f}\frac{\partial^2\mathbf{W}^{(f)}}{\partial t^2}+
\nabla p=-\nabla_{y}\,\Pi
\end{equation}
in the domain $Y_{f}$ for any parameters $(\mathbf{x},t)\in \Omega_{T}$.

There is one equation \eqref{4.18} for two unknown functions
$\mathbf{W}^{(f)}$  and $\Pi$. To derive the second equation we
put in \eqref{4.4}
$\psi=\varepsilon\,h(\mathbf{x},t)\,\psi_0
(\frac{\mathbf{x}}{\varepsilon})$ with arbitrary smooth function $h(\mathbf{x},t)$ and
arbitrary smooth 1-periodic function $\psi_0(\mathbf{y})$ and pass to the limit as
$\varepsilon\to 0$:
\[
\int_{\Omega_{T}}h(\mathbf{x},t)\Big(\int_{Y}\chi(\mathbf{y})
\nabla\psi_0(\mathbf{y})\cdot\mathbf{W}^{(f)}
(\mathbf{x},t,\mathbf{y})dy\Big)dx=0.
\]
After reintegration we obtain the desired microscopic continuity equation
\begin{equation}\label{4.19}
\nabla\cdot\mathbf{W}^{(f)}=0,\quad \mathbf{y}\in Y_{f}.
\end{equation}
A rigorous theory (see \cite{GN1,AM1,AM2}) supplies
the system  \eqref{4.18}, \eqref{4.19} with the boundary condition
\begin{equation}\label{4.20}
\big(\mathbf{W}^{(f)}(\mathbf{x},t,\mathbf{y})-
\mathbf{w}_s(\mathbf{x},t)\big)\cdot\mathbf{n}(\mathbf{y})=0
\end{equation}
on the boundary $\gamma$ with the unit normal
$\mathbf{n}(\mathbf{y})$,  and the homogeneous initial
conditions
\begin{equation}\label{4.21}
\mathbf{W}^{(f)}(\mathbf{x},0,\mathbf{y})=
\frac{\partial\mathbf{W}^{(f)}}{\partial t}(\mathbf{x},0,\mathbf{y})=0.
\end{equation}
Problem \eqref{4.18}-\eqref{4.21} has been solved in \cite{AM2}:
\begin{equation}\label{4.22}
\varrho_{f}\frac{\partial^2\mathbf{W}^{(f)}}{\partial t^2}=
\varrho_{f}\frac{\partial^2\mathbb{w}_s}{\partial t^2}
-\big(\mathbb{I}-\sum_{i=1}^{3}\nabla_{y}\,\Pi_i\otimes\mathbf{e}_i\big)
\cdot\big(\nabla\, p+\varrho_{f}\frac{\partial^2\mathbb{w}_s}{\partial t^2}\big),
\end{equation}
where $\Pi_i(\mathbf{y}),\,i=1,2,3$ are solutions to the periodic boundary value problems
\[
\triangle_{y}\Pi_i=0,\,\,\mathbf{y}\in Y_{f},\,\,
(\nabla_{y}\,\Pi_i-\mathbf{e}_i)\cdot\mathbf{n}(\mathbf{y})=0,
\quad \mathbf{y}\in \gamma.
\]
Thus,
\begin{equation}\label{4.23}
\varrho_{f}\frac{\partial^2\mathbf{w}^{(f)}}{\partial t^2}=
m\varrho_{f}\frac{\partial^2\mathbb{w}_s}{\partial t^2}-\mathbb{B}_2^{(f)}
\cdot\big(\nabla\, p+\varrho_{f}\frac{\partial^2\mathbb{w}_s}{\partial t^2}\big),
\end{equation}
where
\begin{equation}\label{4.24}
\mathbb{B}_2^{(f)}=m\mathbb{I}-\sum_{i=1}^{3}\int_{Y_{f}}
\nabla_{y}\,\Pi_i^{(f)}dy\otimes\mathbf{e}_i.
\end{equation}
Differential equations \eqref{4.6}-\eqref{4.9}, \eqref{4.23},
boundary conditions  \eqref{4.10}-\eqref{4.12}, and initial
conditions \eqref{4.13} form the mathematical model (II) of
seismics in composite media.

\section{Homogenization: case (III)}

According to \cite{AM1} the set of criteria (III) dictates the form
 of the two-scale expansion:
\begin{equation}\label{5.1}
\begin{aligned}
\mathbf{w}^{\varepsilon}(\mathbf{x},t)
&= (1-\varsigma)\mathbf{w}(\mathbf{x},t)+
\varsigma\,\chi(\frac{\mathbf{x}}{\varepsilon})\mathbf{w}^{(f,\varepsilon)}
(\mathbf{x},t)
+\varsigma\big(1-\chi(\frac{\mathbf{x}}{\varepsilon})\big)
\big(\mathbf{w}_s(\mathbf{x},t) \\
&\quad +\varepsilon \mathbf{w}_s^{\varepsilon}(\mathbf{x},t)\big)
+o(\varepsilon),
\\
p^{\varepsilon}(\mathbf{x},t)
&=(1-\varsigma) p(\mathbf{x},t)+
\varsigma\,\chi(\frac{\mathbf{x}}{\varepsilon}) p_{f}(\mathbf{x},t)
+\varsigma\big(1-\chi(\frac{\mathbf{x}}{\varepsilon})\big)
 p^{\varepsilon}_s(\mathbf{x},t)+o(\varepsilon),
\end{aligned}
\end{equation}
where
\[
\mathbf{w}^{(f,\varepsilon)}(\mathbf{x},t)=
\mathbf{W}^{(f)}(\mathbf{x},t,\frac{\mathbf{x}}{\varepsilon}),\,\,\,
\mathbf{w}_s^{\varepsilon}(\mathbf{x},t)=
\mathbf{W}_s(\mathbf{x},t,\frac{\mathbf{x}}{\varepsilon}),\,\,\,
p^{\varepsilon}_s(\mathbf{x},t)=
P_s(\mathbf{x},t,\frac{\mathbf{x}}{\varepsilon}),
\]
and $\mathbf{W}^{(f)}(\mathbf{x},t,\mathbf{y})$,
$\mathbf{W}_s(\mathbf{x},t,\mathbf{y})$,
$P_s(\mathbf{x},t,\mathbf{y})$ are 1-periodic in the variable
 $\mathbf{y}$ functions.

Next we express the pressure $p^{\varepsilon}_s$ in the solid
component  in $\Omega$ using the continuity equation \eqref{2.8} and
two-scale expansion \eqref{5.1}:
\begin{equation}\label{5.2}
p^{\varepsilon}_s(\mathbf{x},t)
=-\bar{c}_s^2\big( \nabla\cdot \mathbf{w}_s(\mathbf{x},t)+\nabla_{y}\cdot
\mathbf{W}_s(\mathbf{x},t,\frac{\mathbf{x}}{\varepsilon})\big)+o(\varepsilon).
\end{equation}
The substitution of \eqref{5.1} and \eqref{5.2} into \eqref{2.10} results 
in the integral equality
\begin{equation}\label{5.3}
\begin{aligned}
&\int_{Q_{T}}\nabla\cdot(p^0\mathbf{\varphi})\,dx\,dt+
\int_{\Omega^{(0)}_{T}}\big(\varrho^{(0)}_s
\frac{\partial^2\mathbf{w}}{\partial t^2}\cdot
\mathbf{\varphi}+\big(\lambda_0^{(0)}\mathbb{D}(x,\mathbf{w})-
p\mathbb{I}\big):\mathbb{D}(x,\mathbf{\varphi})\big)\,dx\,dt\\
&+\int_{\Omega_{T}}\Big(\big(\varrho_{f}\chi(\frac{\mathbf{x}}{\varepsilon})
\frac{\partial^2\mathbf{W}^{(f)}}{\partial t^2}
(\mathbf{x},t,\frac{\mathbf{x}}{\varepsilon})
+\varrho_s\big(1-\chi(\frac{\mathbf{x}}{\varepsilon})\big)
\frac{\partial^2\mathbf{w}_s}{\partial t^2}\big)
\cdot\mathbf{\varphi}\Big)\,dx\,dt\\
&+\int_{\Omega_{T}}\lambda_0\big(1-\chi(\frac{\mathbf{x}}{\varepsilon})\big)
\Big(\mathfrak{N}^{(0)}:\big(\mathbb{D}(x,\mathbf{w}_s)+\mathbb{D}
\big(y,\mathbf{W}_s(\mathbf{x},t,\frac{\mathbf{x}}{\varepsilon})\big)
\big)\Big):\mathbb{D}(x,\mathbf{\varphi})\,dx\,dt\\
&= -\int_{\Omega_{T}}\chi(\frac{\mathbf{x}}{\varepsilon})\Big(\frac{\alpha_{\mu}}
{\varepsilon}\mathbb{D}\big(y,\frac{\partial\mathbf{W}^{(f)}}{\partial t}
(\mathbf{x},t,\frac{\mathbf{x}}{\varepsilon})\big)
- p_{f}\mathbb{I}\Big):\mathbb{D}(x,\mathbf{\varphi})\,dx\,dt+o(\varepsilon),
\end{aligned}
\end{equation}
which holds for any smooth function $\mathbf{\varphi}(\mathbf{x},t)$.
In \eqref{5.3}
\[
\mathfrak{N}^{(0)}=\sum_{i,j=1}^{3}\mathbb{J}^{ij}\otimes\mathbb{J}^{ij}+
\frac{\bar{c}_s^2}{\lambda_0}\mathbb{I}\otimes\mathbb{I},\,\,\,
\mathbb{J}^{ij}=\frac{1}{2}\big(\mathbf{e}_i\otimes\mathbf{e}_{j}+
\mathbf{e}_{j}\otimes\mathbf{e}_i\big),
\]
$\{\mathbf{e}_1,\mathbf{e}_2,\mathbf{e}_3\}$ is a
standard  Cartesian basis, and the fourth-rank tensor
$\mathbb{A}\otimes\mathbb{B}$ is the tensor (direct) product of the
second-rank tensors $\mathbb{A}$ and $\mathbb{B}$:
\[
(\mathbb{A}\otimes\mathbb{B}):\mathbb{C}=\mathbb{A}(\mathbb{B}:\mathbb{C})
\]
for any second-rank tensor $\mathbb{C}$.

After the limit in \eqref{5.3} as $\varepsilon \to 0$ we arrive at the 
integral identity
\begin{equation}\label{5.4}
\begin{aligned}
&\int_{\Omega^{(0)}_{T}}\big(\varrho^{(0)}_s
\frac{\partial^2\mathbf{w}}{\partial t^2}\cdot\mathbf{\varphi}
+\big(\lambda_0^{(0)}\mathbb{D}(x,\mathbf{w})-p\mathbb{I}\big):
\mathbb{D}(x,\mathbf{\varphi})\big)\,dx\,dt\\
&+\int_{Q_{T}}\nabla\cdot(p^0\mathbf{\varphi})\,dx\,dt
+\int_{\Omega_{T}}\Big(\big(\varrho_{f}
\frac{\partial^2\mathbf{w}^{(f)}}{\partial t^2}
+\varrho_s(1-m)\frac{\partial^2\mathbf{w}_s}{\partial t^2}\big)
\cdot\mathbf{\varphi}\Big)\,dx\,dt\\
&+\int_{\Omega_{T}}\lambda_0\Big(\mathfrak{N}^{(0)}:\big((1-m)
\mathbb{D}(x,\mathbf{w}_s)+
\langle\mathbb{D}(y,\mathbf{W}_s)\rangle_{Y_s}\big)\Big):
\mathbb{D}(x,\mathbf{\varphi})\,dx\,dt\\
&=\int_{\Omega_{T}}\big(m p_{f}\mathbb{I}\big):\mathbb{D}(x,\mathbf{\varphi})
\,dx\,dt,
\end{aligned}
\end{equation}
where
\[
\mathbf{w}^{(f)}=\langle\,\mathbf{W}^{(f)}\,\rangle_{Y_{f}},\quad
\langle F\rangle_{A}
=\int_{A}\,F(\mathbf{y})dy,\quad A\subseteq Y.
\]
To pass to the limit as $\varepsilon \to 0$ in the
continuity equation  \eqref{2.8} we rewrite it as an integral
identity and use the representation \eqref{5.1}:
\begin{equation}\label{5.5}
\begin{aligned}
&\int_{Q_{T}}\psi\,\Big(\frac{(1-\varsigma)}{\bar{c}_{s,0}^2} p+
\varsigma\,\chi(\frac{\mathbf{x}}{\varepsilon})\frac{p_{f}}{\bar{c}_{f}^2}+
\varsigma\big(1-\chi(\frac{\mathbf{x}}{\varepsilon})
\frac{p_s^{\varepsilon}}{\bar{c}_s^2}\big)\Big)\,dx\,dt\\
&-\int_{Q_{T}}\nabla\,\psi\cdot\big((1-\varsigma)\mathbf{w}+
\varsigma\,\chi(\frac{\mathbf{x}}{\varepsilon})\mathbf{W}_{f}+
\varsigma\big(1-\chi(\frac{\mathbf{x}}{\varepsilon})\big)
\mathbf{w}_s\big)\,dx\,dt=o(\varepsilon).
\end{aligned}
\end{equation}
In the limit as $\varepsilon \to 0$  results in the integral equality
\begin{equation}\label{5.6}
\begin{aligned}
&\int_{Q_{T}}\psi\big(\frac{(1-\varsigma)}{\bar{c}_{s,0}^2} p+
\frac{\varsigma\,m}{\bar{c}_{f}^2} p_{f}+
\frac{\varsigma}{\bar{c}_s^2}\,\langle\,P_s
\rangle_{Y_s}\big))\,dx\,dt\\
&-\int_{Q_{T}}\nabla\,\psi\cdot\big((1-\varsigma)\mathbf{w}+
\varsigma\,\mathbf{w}^{(f)}+\varsigma(1-m)\mathbf{w}_s\big)\,dx\,dt=0,
\end{aligned}
\end{equation}
which holds for any smooth function $\psi$ vanishing at $S$.

Finally, we rewrite the continuity equation in the pore space
 $\Omega_{f}$ as the corresponding integral identity
\begin{align*}
&\int_{\Omega_{T}}\psi\big(\frac{\chi^{\varepsilon}}{\bar{c}_{f}^2}
p^{\varepsilon}+
\chi^{\varepsilon}\nabla\cdot\mathbf{w}^{\varepsilon}\big)\,dx\,dt\\
&= \int_{\Omega_{T}}\psi\big(\frac{\chi^{\varepsilon}}{\bar{c}_{f}^2}
 p^{\varepsilon}+ \nabla\cdot\mathbf{w}^{\varepsilon}-(1-\chi^{\varepsilon})\,
\nabla\cdot\mathbf{w}^{\varepsilon}\big)\,dx\,dt\\
&=\int_{\Omega_{T}}\big(\psi\,\frac{\chi^{\varepsilon}}{\bar{c}_{f}^2}
 p^{\varepsilon}- (\nabla\psi)\cdot\mathbf{w}^{\varepsilon}
 -\psi(1-\chi^{\varepsilon})
\nabla\cdot\mathbf{w}^{\varepsilon}\big)\,dx\,dt,
\end{align*}
and apply the two-scale expansion \eqref{5.1}:
\begin{equation}\label{5.7}
\begin{aligned}
&\int_{\Omega_{T}}\Big(\frac{\psi}{\bar{c}_{f}^2}\,
\chi(\frac{\mathbf{x}}{\varepsilon}) p_{f}-(\nabla\psi)\cdot
\big(\chi(\frac{\mathbf{x}}{\varepsilon})\mathbf{W}_{f}^{\varepsilon}+
\big(1-\chi(\frac{\mathbf{x}}{\varepsilon})\big)\mathbf{w}_s\big)\\
&-\psi\big(1-\chi(\frac{\mathbf{x}}{\varepsilon})\big)
\big(\nabla\cdot\mathbf{w}_s+
\nabla_{y}\cdot\mathbf{W}_s\big)\Big)\,dx\,dt=o(\varepsilon).
\end{aligned}
\end{equation}
In the limit as $\varepsilon \to 0$ results in the desired integral equality
\begin{equation}\label{5.8}
\int_{\Omega_{T}}\big(\psi\,\frac{m}{\bar{c}_{f}^2} p_{f}-
\nabla\,\psi\cdot\mathbf{w}^{(f)}-
\psi\,\langle\nabla_{y}\cdot\mathbf{W}_s\rangle_{Y_s}\big)\,dx\,dt=0.
\end{equation}
The localization of \eqref{5.4}, \eqref{5.6}, and \eqref{5.8} gives as 
the Lame system
\begin{gather}\label{5.9}
\varrho^{(0)}_s\frac{\partial^2\mathbf{w}}{\partial t^2}=
\nabla\cdot \mathbb{P}^{(0)}_s,\quad
\mathbb{P}^{(0)}_s=\lambda^{(0)}_0\,
\mathbb{D}(x,\mathbf{w})-p\mathbb{I}, \\
\label{5.10}
p+\bar{c}_{s,0}^2\nabla\cdot \mathbf{w}=0
\end{gather}
in the domain $\Omega^{(0)}$ for $t>0$, the macroscopic dynamic equation
\begin{gather}\label{5.11}
\varrho_{f}\frac{\partial^2\mathbf{w}^{(f)}}{\partial t^2}+
\varrho_s(1-m)\frac{\partial^2\mathbf{w}_s}{\partial t^2}=
\nabla\cdot \widehat{\mathbb{P}}, \\
\label{5.12}
\widehat{\mathbb{P}}=\lambda_0\,\mathfrak{N}^{(0)}:
\big((1-m)\mathbb{D}(x,\mathbf{w}_s)+
\langle\mathbb{D}(y,\mathbf{W}_s)\rangle_{Y_s}\big)-m p_{f}\mathbb{I}
\end{gather}
for the solid component and the macroscopic continuity equation
\begin{equation}\label{5.13}
\frac{m}{\bar{c}_{f}^2} p_{f}+\nabla\cdot\mathbf{w}^{(f)}=
\langle\nabla_{y}\cdot\mathbf{W}_s\rangle_{Y_s}
\end{equation}
for the liquid component in the domain $\Omega$ for $t>0$.

The same localization of \eqref{5.4} and \eqref{5.6} also provides 
the boundary condition
\begin{equation}\label{5.14}
\mathbb{P}^{(0)}_s\cdot\mathbf{e}_3=p^0\cdot\mathbf{e}_3
\end{equation}
on the outer boundary $S$ with the unit normal $\mathbf{e}_3$, 
and the continuity conditions
\begin{gather}\label{5.15}
\lim_{\substack{\mathbf{x}\to \mathbf{x}^0\\
\mathbf{x}\in \Omega^{(0)}}}\mathbb{P}^{(0)}(\mathbf{x},t)
\cdot\mathbf{e}_3=
\lim_{\substack{\mathbf{x}\to \mathbf{x}^0\\
\mathbf{x}\in \Omega}}\widehat{\mathbb{P}}(\mathbf{x},t)
\cdot\mathbf{e}_3,
\\
\label{5.16}
\lim_{\substack{\mathbf{x}\to \mathbf{x}^0\\
\mathbf{x}\in \Omega^{(0)}}}
\mathbf{w}(\mathbf{x},t)\cdot\mathbf{e}_3=
\lim_{\substack{\mathbf{x}\to \mathbf{x}^0\\
\mathbf{x}\in \Omega}}\big(\mathbf{w}^{(f)}(\mathbf{x},t)+
(1-m)\mathbf{w}_s(\mathbf{x},t)\big)\cdot\mathbf{e}_3
\end{gather}
on the common boundary $S^{(0)}\ni\mathbf{x}^0$ with the unit normal
$\mathbf{e}_3$.

More detailed mathematical analysis shows that
\begin{equation}\label{5.17}
\lim_{\substack{\mathbf{x}\to \mathbf{x}^0\\
\mathbf{x}\in \Omega^{(0)}}}
\mathbf{w}(\mathbf{x},t)=
\lim_{\substack{\mathbf{x}\to \mathbf{x}^0\\
\mathbf{x}\in \Omega}}
(1-m)\mathbf{w}_s(\mathbf{x},t)
\end{equation}
for $\mathbf{x}^0\in S^{(0)}$. Unfortunately we have no
possibility  to prove  the statement due  technical reasons.

Differential equations and boundary conditions are supplemented with 
initial conditions
\begin{equation}\label{5.18}
\begin{aligned}
\mathbf{w}(\mathbf{x},0)
&=\frac{\partial\mathbf{w}}{\partial t}(\mathbf{x},0)=
\mathbf{w}_s(\mathbf{x},0)=
\frac{\partial\mathbf{w}_s}{\partial t}(\mathbf{x},0)\\
&=\mathbf{w}^{(f)}(\mathbf{x},0)=
\frac{\partial\mathbf{w}^{(f)}}{\partial t}(\mathbf{x},0)=0.
\end{aligned}
\end{equation}
However, the obtained system \eqref{5.9}-\eqref{5.18} is still
incomplete. We need two more  differential equations for
$\mathbf{W}_s$ and $\mathbf{W}^{(f)}$. More precisely, we
have to express the terms
$\langle\mathbb{D}(y,\mathbf{W}_s)\rangle_{Y_s}$ and
$\langle\nabla_{y}\cdot\mathbf{W}_s\rangle_{Y_s}$ by means
of functions $\mathbb{D}(x,\mathbf{w}_s)$ and $p_{f}$ and
rewrite \eqref{5.12} and \eqref{5.13} as
\begin{gather}\label{5.19}
\widehat{\mathbb{P}}=\lambda_0\,\mathfrak{N}^{s}:\mathbb{D}(x,\mathbf{w}_s)-
 p_{f}\mathbb{C}^{s}, \\
\label{5.20}
\frac{m}{\bar{c}_{f}^2} p_{f}+\nabla\cdot\mathbf{w}^{(f)}=
\mathbb{C}_0^{s}:\mathbb{D}(x,\mathbb{w}_s)+\frac{c^{s}_0}{\lambda _0}p_{f}.
\end{gather}
To find the missing equation for the function $\mathbf{W}_s$ let
us consider the integral  identity \eqref{5.3}.  As in previous
section, we choose test functions in the form
$\mathbf{\varphi}^{\varepsilon}
(\mathbf{x},t)=\varepsilon\,h(\mathbf{x},t)\,\mathbf{\varphi}_0
(\frac{\mathbf{x}}{\varepsilon})$, where $h$ is an arbitrary
smooth function with a compact support in $\Omega$ vanishing on $S$,
and $\mathbf{\varphi}_0(\mathbf{y})$ is an arbitrary
1-periodic smooth function with a compact support in $Y_s$.

The limit in \eqref{5.3} as $\varepsilon \to 0$ with chosen test functions results in
\begin{align*}
&\int_{\Omega_{T}}h\,\Big(\int_{Y}\Big(\lambda_0\big(1-\chi(\mathbf{y})\big)
\big(\mathfrak{N}^{(0)}:(\mathbb{D}(x,\mathbf{w}_s)\\
&+\mathbb{D}(y,\mathbf{W}_s)\big)
-\chi(\mathbf{y})\,m p_{f}\mathbb{I}\Big):
\mathbb{D}(y,\mathbf{\varphi}_0)dy\Big)\,dx\,dt=0
\end{align*}
After a localization we obtain the differential equation
\begin{equation}\label{5.21}
\nabla_{y}\cdot\Big(\lambda_0\big(1-\chi(\mathbf{y})\big)
\mathfrak{N}^{(0)}:\big(\mathbb{D}(x,\mathbf{w}_s)+
\mathbb{D}(y,\mathbf{W}_s)\big)-m p_{f}\,\chi(\mathbf{y})\Big)=0
\end{equation}
in the domain $Y$, which is understood in the sense of
distributions.  That is, as a usual differential equation
\begin{equation}\label{5.22}
\nabla_{y}\cdot\Big(\mathfrak{N}^{(0)}:\big(\mathbb{D}(x,\mathbf{w}_s)+
\mathbb{D}(y,\mathbf{W}_s)\big)\Big)=0
\end{equation}
in the domain $Y_s$. In the same way using test functions with a compact
support localizes at $\gamma$ we derive the boundary condition
\begin{equation}\label{5.23}
\Big(\lambda_0\mathfrak{N}^{(0)}:\big(\mathbb{D}(x,\mathbf{w}_s)+
\mathbb{D}(y,\mathbf{W}_s)\big)\Big)\cdot\mathbf{n}=
-m p_{f}\,\mathbf{n}
\end{equation}
on the boundary $\gamma$. Here $\mathbf{n}$ is a unit normal to $\gamma$.

The problem \eqref{5.18}, \eqref{5.19} is completed with  the
periodicity conditions on the remaining part $\partial
Y_s\backslash\gamma$ of the boundary $\partial Y_s$.

Let $\mathbf{U}_2^{(ij)}(\mathbf{y})$ and $\mathbf{U}_2^{(0)}(\mathbf{y})$
be solutions of periodic problems
\begin{gather} \label{5.24}
\nabla_{y}\cdot\Big((1-\chi)\Big(\mathfrak{N}^{(0)}:\big(\mathbb{J}^{(ij)}+
\mathbb{D}(y,\mathbf{U}_2^{(ij)})\big)\Big)\Big)=0, \\
\label{5.25}
\nabla_{y}\cdot\Big((1-\chi)\big(\mathfrak{N}^{(0)}:
\mathbb{D}(y,\mathbf{U}_2^{(0)})+\mathbb{I}\big)\Big)=0
\end{gather}
in $Y$.
Then
\[
\mathbf{W}_s(\mathbf{x},t,\mathbf{y})=
\sum_{i,j=1}^{3}\mathbf{U}_2^{(ij)}(\mathbf{y})D_{ij}(\mathbf{x},t)+
\frac{m}{\lambda_0} p_{f}(\mathbf{x},t)\,\mathbf{U}_2^{(0)}(\mathbf{y}),
\]
where
\[
D_{ij}=\frac{1}{2}\Big(\frac{\partial u_i}{\partial x_{j}}+
\frac{\partial u_{j}}{\partial x_i}\Big),
\quad \mathbb{w}_s=(u_1,\,u_2,\,u_3),\quad
\mathbb{D}(x,\mathbf{w}_s)=\sum_{i,j=1}^{3}D_{ij}\mathbb{J}^{(ij)}.
\]
Thus
\begin{align*}
&\langle\mathbb{D}(y,\mathbf{W}_s)\rangle_{Y_s}\\
&=\sum_{i,j=1}^{3}\langle\mathbb{D}(y,\mathbf{U}_2^{(ij)})\rangle_{Y_s}D_{ij}+
\frac{m}{\lambda_0} p_{f}
\langle\mathbb{D}(y,\mathbf{U}_2^{(0)})\rangle_{Y_s}\\
&=\Big(\sum_{i,j=1}^{3}\langle\mathbb{D}(y,\mathbf{U}_2^{(ij)})\rangle_{Y_s}
\otimes\mathbb{J}^{(ij)}\Big):\mathbb{D}(x,\mathbb{w}_s)+\frac{m}{\lambda_0}
 p_{f}\,\langle\mathbb{D}(y,\mathbf{U}_2^{(0)})\rangle_{Y_s},
\end{align*}
\begin{align*}
&\lambda _0\,\mathfrak{N}^{(0)}:\big((1-m)\mathbb{D}(x,\mathbb{w}_s)+
\langle\mathbb{D}(y,\mathbf{W}_s)\rangle_{Y_s}\big)-m p_{f}\mathbb{I}\\
&=\lambda _0\,\mathfrak{N}^{s}:\mathbb{D}(x,\mathbb{w}_s)-
p_{f}\mathbb{C}^{s},
\end{align*}
\begin{align*}
\langle\nabla_{y}\cdot\mathbf{W}_s\rangle_{Y_s}
&= \sum_{i,j=1}^{3}\langle\nabla_{y}\cdot\mathbf{U}_2^{(ij)}\rangle_{Y_s}D_{ij}+
\frac{m}{\lambda_0} p_{f}\langle\nabla_{y}\cdot
\mathbf{U}_2^{(0)}\rangle_{Y_s} \\
&= \big(\sum_{i,j=1}^{3}\langle\nabla_{y}\cdot
\mathbf{U}_2^{(ij)}\rangle_{Y_s}\mathbb{J}^{ij}\big):
\mathbb{D}(x,\mathbb{w}_s)+\big(\frac{m}{\lambda_0}\langle\nabla_{y}
\cdot\mathbf{U}_2^{(0)}\rangle_{Y_s}\big) p_{f},
\end{align*}
where
\begin{gather}\label{5.26}
\mathfrak{N}^{s}=\mathfrak{N}^{(0)}:
\Big((1-m)\sum_{i,j=1}^{3}\mathbb{J}^{ij}\otimes\mathbb{J}^{ij}+
\sum_{i,j=1}^{3}\langle\mathbb{D}(y,\mathbf{U}_2^{(ij)})\rangle_{Y_s}
\otimes\mathbb{J}^{(ij)}\Big), \\
 \label{5.27}
\mathbb{C}^{s}=m\mathbb{I}-\langle\mathbb{D}(y,\mathbf{U}_2^{(0)})\rangle_{Y_s}, \\
 \label{5.28}
\mathbb{C}^{s}_0=\sum_{i,j=1}^{3}\langle\nabla_{y}
\cdot\mathbf{U}_2^{(ij)}\rangle_{Y_s}\mathbb{J}^{ij},\,\,
c^{s}_0=\langle\nabla_{y}\cdot\mathbf{U}_2^{(0)}\rangle_{Y_s}.
\end{gather}
The derivation of the equation for $\mathbf{W}^{(f)}$ repeats in
its main  features the arguments of the previous section. We choose the
test functions $\mathbf{\varphi}^{\varepsilon}$ in \eqref{5.3}
as
\[
\mathbf{\varphi}^{\varepsilon}(\mathbf{x},t)=h(\mathbf{x},t)\,
\mathbf{\varphi}_0(\frac{\mathbf{x}}{\varepsilon}),
\]
where $h$ is a smooth function with a compact support in $\Omega$
and $\mathbf{\varphi}_0(\mathbf{y})$ is a smooth
1-periodic solenoidal function with a compact support in $Y_{f}$.
After the limit as $\varepsilon\to 0$ and localisation we
arrive at the differential equation
\begin{equation}\label{5.29}
\varrho_{f}\frac{\partial^2\mathbf{W}^{(f)}}{\partial t^2}=
\mu_1\nabla\cdot\mathbb{D}\big(y,\frac{\partial\mathbf{W}^{(f)}}{\partial t}\big)-
\nabla_{y}\Pi^{(f)}-\nabla p_{f}
\end{equation}
in the domain $Y_{f}$ for $t>0$. Here, as in previous section, we
also  must define a 1-periodic in $\mathbf{y}$ function
$\Pi^{(f)}(\mathbf{x},t,\mathbf{y})$, which appears due to
the choice of test functions. The missing equation is derived from
the continuity equation in its integral form \eqref{5.5} in the same
way as in the previous section and coincides with \eqref{4.19}.

According to \cite{AM1} and \cite{AM2} the system \eqref{5.29}, \eqref{4.19} supplies  with the boundary condition
\begin{equation}\label{5.30}
\mathbf{W}^{(f)}(\mathbf{x},t,\mathbf{y})=
\mathbf{w}_s(\mathbf{x},t)
\end{equation}
on the boundary $\gamma$, and the homogeneous initial conditions
\begin{equation}\label{5.31}
\mathbf{W}^{(f)}(\mathbf{x},0,\mathbf{y})=
\frac{\partial\mathbf{W}^{(f)}}{\partial t}(\mathbf{x},0,\mathbf{y})=0.
\end{equation}
Problem \eqref{4.19}, \eqref{5.29}-\eqref{5.31} has been solved in \cite{AM2}:
\begin{align*}
\mathbf{W}^{(f)}
&=\mathbb{w}_s(\mathbf{x},t)+
\sum_{i=1}^{3}\int_0^{t}\mathbf{W}^{(f)}_i(\mathbf{y},t-\tau)
\big(\frac{\partial p_{f}}{\partial x_i}(\mathbf{x},\tau)+
\varrho_{f}\frac{\partial^2{w}_{s,i}}{\partial \tau^2}
(\mathbf{x},\tau)\big)d\tau \\
&=\mathbb{w}_s(\mathbf{x},t)+
\sum_{i=1}^{3}\int_0^{t}\big(\mathbf{W}^{(f)}_i(\mathbf{y},t-\tau)
\otimes\mathbf{e}_i\big)\cdot\big(\nabla p_{f}(\mathbf{x},\tau)+
\varrho_{f}\frac{\partial^2\mathbf{w}_s}{\partial \tau^2}
(\mathbf{x},\tau)\big)d\tau,
\end{align*}
\[
\Pi^{(f)}(\mathbf{x},t,\mathbf{y})=\sum_{i=1}^{3}\int_0^{t}\Pi^{(f)}_i
(\mathbf{y},t-\tau)\big(\frac{\partial p_{f}}{\partial x_i}(\mathbf{x},\tau)+
\varrho_{f}\frac{\partial^2{w}_{s,i}}{\partial \tau^2}
(\mathbf{x},\tau)\big)d\tau,
\]
where $\mathbb{w}_s=({w}_{s,1},{w}_{s,2},{w}_{s,3})$ and
$\{\mathbf{W}^{(f)}_i,\,\Pi^{(f)}_i\}$ , $i=1,2,3$,  are
solutions to the following periodic initial boundary value problem
\begin{gather}\label{5.32}
\varrho_{f}\frac{\partial^2\mathbf{W}^{(f)}_i}{\partial t^2}
=\mu_1\nabla\cdot\mathbb{D}\big(y,
\frac{\partial\mathbf{W}^{(f)}_i}{\partial t})-
\nabla_{y}\Pi^{(f)}_i,\quad(\mathbf{y},t)\in Y_{f}\times(0,T), \\
\label{5.33}
\nabla_{y}\cdot\mathbf{W}^{(f)}_i(\mathbf{y},t)=0,
\,(\mathbf{y},t)\in Y_{f}\times(0,T), \\
\label{5.34}
\mathbf{W}^{(f)}_i(\mathbf{y},0)=0,\,\,
\varrho_{f}\frac{\partial\mathbf{W}^{(f)}_i}{\partial t}(\mathbf{y},0)
=-\mathbf{e}_i,\quad \mathbf{y}\in Y_{f}, \\
\label{5.35}
\mathbf{W}^{(f)}_i(\mathbf{y},t)=0,\,(\mathbf{y},t)\in \gamma\times(0,T).
\end{gather}
Thus,
\begin{equation}\label{5.36}
\begin{aligned}
\mathbf{w}^{(f)}(\mathbf{x},t)
&=\int_{Y_{f}}\mathbf{W}^{(f)}
(\mathbf{x},t,\mathbf{y})dy\\
&=m\,\mathbb{w}_s(\mathbf{x},t)
+\int_0^{t}\mathbb{B}_3^{(f)}(t-\tau)\cdot\big(\nabla p(\mathbf{x},\tau)+
\varrho_{f}\frac{\partial^2\mathbf{w}_s}{\partial \tau^2}
(\mathbf{x},\tau)\big)d\tau,
\end{aligned}
\end{equation}
where
\begin{equation}\label{5.37}
\mathbb{B}_3^{(f)}(t)=\sum_{i=1}^{3}\int_{Y_{f}}\mathbf{W}^{(f)}_i
(\mathbf{y},t)dy\otimes\mathbf{e}_i.
\end{equation}
Differential equations \eqref{5.9}-\eqref{5.11}, \eqref{5.20}, and
\eqref{5.36},  boundary conditions \eqref{5.14}-\eqref{5.17},
initial conditions \eqref{5.18} and state equations  \eqref{5.19}
and  \eqref{5.26}-\eqref{5.28} constitute the mathematical
model (III) of seismics in composite media.

\section{One dimensional model for the case (I): numerical implementations}

\subsection*{Direct problem}
For the sake of simplicity we consider the space,  which consists of
the following subdomains: 
$\Omega_1=\{x\in \mathbb{R}:0<x<H_1\}$,
$\Omega_2=\{x\in \mathbb{R}: H_1<x<H_2\}$, and
$\Omega_3=\{x\in \mathbb{R}:x>H_2\}$. Differential equations
\eqref{3.8}, \eqref{3.9}, \eqref{3.12}, and \eqref{3.13}
result in
\[
\frac{1}{\hat{c}^2(x)}\frac{\partial^2 p}{\partial
t^2}=\operatorname{div}\big(\frac{1}{\hat{\rho}(x)}\nabla
(p+m\frac{\nu_0}{\bar{c}_{f}^2}\frac{\partial p}{\partial t})\big)
\]
where
\[
\frac{1}{\hat{c}^2}=\frac{m}{\bar{c}_{f}^2}+\frac{(1-m)}{\bar{c}_s^2}
\]
and $\hat{\rho}=m\rho_{f}+(1-m)\rho_s$ are respectively average
wave propagation velocity and average density of the medium.

Applying now the Fourier transformation we arrive at
\begin{equation}\label{6.1}
\frac{d^2\hat{P}}{dX^2}+\frac{\hat{\rho}\omega^2}
{(1-\frac{m\nu_0}{c^2_{f}}i\omega)\hat{c}^2}\hat{P}=0
\end{equation}
where
$\hat{P}{(x, \omega)}$-the pressure obtained after Fourier transform.

\begin{figure}[ht]
\begin{center}
 \includegraphics[width=0.7\textwidth]{fig3} maratainur.jpg
\end{center}
\caption{Scheme of arrangement of layers}
 \label{fig3}
\end{figure}

Depending on the exact physical properties, the sedimentary rock
zone is divided into three subdomains. The value of the geometry of
pores, viscosity of fluid, density of rock, and velocity of seismic
wave considered in each layers to be different. In the experiment in
order to get numerical solution, it's assumed that the first medium
is a shale, the second medium is oil saturated sandstone, and the
third medium is a limestone (see Fig.\ref{fig3}).

Let us suppose that there is a plane wave which propagates from
$\infty$. Then the general solution of equation \eqref{6.1} for
$-\infty<X\leq H_1$ in the case $\nu_0=0$ is written down as:
\begin{equation}\label{6.2}
\hat{P_1}=\exp\Big\{\frac{i\omega\sqrt{\hat{\rho_1}}}{\hat{c_1}}x\Big\}
+A_2\exp\Big\{\frac{-i\omega\sqrt{\hat{\rho_1}}}{\hat{c_1}}x\Big\}.
\end{equation}
The general solution of equation \eqref{6.1} for $H_1\leq x<H_2$ in
the case $\nu_0>0$ is represented as:
\begin{equation}\label{6.3}
\hat{P_2}=B_1\exp\Big\{\frac{i\omega\sqrt{\hat{\rho_2}}}{\hat{c}_2
\sqrt{1-\frac{m\nu_0}{c^2_{f}}i\omega}}x\big\}
+B_2\exp\Big\{\frac{-i\omega\sqrt{\hat{\rho_2}}}{\hat{c}_2
\sqrt{1-\frac{m\nu_0}{c^2_{f}}i\omega}}x\Big\}.
\end{equation}
Finally the general solution for $x\geq H_2$ in the case $\nu_0=0$
will be the following:
\begin{equation}\label{6.4}
\hat{P}_3=D_1\exp\Big\{i\omega\frac{\sqrt{\hat{\rho}_3}}{\hat{c}_3}x\Big\}.
\end{equation}
Continuity condition in contact media will be:
\begin{gather}\label{6.5}
[\hat{P_1}-i\omega\frac{m\nu_0}{c^2_{f}}\hat{P_1}]_{H_{1-0}}
=[\hat{P_2}-i\omega\frac{m\nu_0}{c^2_{f}}\hat{P_2}]_{H_{1+0}} \\
\label{6.6}
\hat{c}^2_1\frac{d}{dx}[\hat{P_1}
-i\omega\frac{m\nu_0}{c^2_{f}}\hat{P_1}]_{H_{1-0}}
=\hat{c}^2_2\frac{d}{dx}[\hat{P_2}-i\omega\frac{m\nu_0}{c^2_{f}}
\hat{P_2}]_{H_{1+0}} \\
\label{6.7}
[\hat{P_2}-i\omega\frac{m\nu_0}{c^2_{f}}\hat{P_2}]_{H_{2-0}}
=[\hat{P_3}-i\omega\frac{m\nu_0}{c^2_{f}}\hat{P_3}]_{H_{2+0}} \\
\label{6.8}
\hat{c}^2_2\frac{d}{dx}[\hat{P_2}
-i\omega\frac{m\nu_0}{c^2_{f}}\hat{P_2}]_{H_{2-0}}
=\hat{c}^2_3\frac{d}{dx}[\hat{P_3}
-i\omega\frac{m\nu_0}{c^2_{f}}\hat{P_3}]_{H_{2+0}}
\end{gather}
These relations are nothing else but the system of linear algebraic
equations for the coefficients $A_2, B_1, B_2, D_1$ which
can be easily resolved by any direct method. These coefficients are
used in order to construct the solution in time frequency domain and
after inverse Fourier transform in time the solution in the time
domain can be easily recovered (see Fig.\ref{fig4}).

\begin{figure}[ht]
\begin{center}
 \includegraphics[width=0.7\textwidth]{fig4} % direct.jpg
\end{center}
\caption{Propagation of seismic waves in different layers}
\label{fig4}
\end{figure}

\subsection*{Inverse problem}

In inverse problem \cite{A-F} except $\hat{P}(x,\omega)$ the values
$H_1$, $H_2$, $\hat{c}_2$, $\nu_0$, $m$ are unknown as well.
To determine these values one needs some additional information
about solution of the direct problem - data of inverse problem.
Usually they are given as function $\overline{P}(\omega)$ at $X=0$.
The most widespread way is to search for these values by
minimization of the data misfit functional being $L_2$ norm of the
difference of given functions and computed for some current values
of unknown parameters:
\begin{equation}\label{6.9}
F_i(H^i_1, H^i_2)=\int_{\omega_1}^{\omega_n}
|\hat{P}_i(\omega, H^i_1, H^i_2)-\overline{P}(\omega,
H_1, H_2)|^2 d\omega\to 0
\end{equation}
\begin{equation}\label{6.10}
F_i(H^i_1, \hat{c}^i_2)=\int_{\omega_1}^{\omega_n}
|\hat{P}_i(\omega, H^i_1,
\hat{c}^i_2)-\overline{P}(\omega, H_1, \hat{c}_2)|^2
d\omega\to 0
\end{equation}
\begin{equation}\label{6.11}
F_i(H^i_2, \hat{c}^i_2)=\int_{\omega_1}^{\omega_n}
|\hat{P}_i(\omega, H^i_2,
\hat{c}^i_2)-\overline{P}(\omega, H_2, \hat{c}_2)|^2
d\omega\to 0
\end{equation}
Here $\overline{P}(\omega,\ldots,\ldots)$ is the given wave fields at
$X=0$, while $\hat{P}(\omega,\ldots,\ldots)$ are wave fields
computed for some current values of the desired parameters.

In our numerical experiments the minimum is searched by the
Nelder-Mead technique  (\cite{M-F}, Fig.\ref{fig5}).

\begin{figure}[htbp]
\begin{center}
 \includegraphics[width=0.4\textwidth]{fig5} %neldermethod.jpg
\end{center}
\caption{Simple scheme of Nelder-Mead for two variables regular simplex }
\label{fig5}
\end{figure}

\subsubsection*{Recovery of $H_1$ and $H_2$}
Behavior of the data misfit functional for this statement is
represented in Figures \ref{fig6} and \ref{fig7}. As one
can see this functional is convex and has the unique minimum point.
Therefore this inverse problem is well resolved.

\begin{figure}[ht]
\begin{center}
 \includegraphics[width=0.7\textwidth]{fig6} surfWsveteng.jpg
\end{center}
\caption{Minimization of the functional $F(H_1, H_2)$.}
\label{fig6}
\end{figure}

\begin{figure}[htbp]
\begin{center}
 \includegraphics[width=0.7\textwidth]{fig7} contour200eng.jpg
\end{center}
\caption{Level line of the functional $F(H_1,H_2)$.}
\label{fig7}
\end{figure}

\subsubsection*{Recovery of $H_1$ and $c_2$}
Now we come to the non convex functional and therefore inverse
problem may have few solutions (see Figures \ref{fig8} and \ref{fig9}).

\begin{figure}[ht]
\begin{center}
 \includegraphics[width=0.7\textwidth]{fig8} % surfW01001sveteng.jpg
\end{center}
\caption{Minimization of the functional $F(H_1, \hat{c}_2)$.}
\label{fig8}
\end{figure}


\begin{figure}[ht]
\begin{center}
 \includegraphics[width=0.7\textwidth]{fig9} % 50contour01001eng.jpg
\end{center}
\caption{Level line functional $F(H_1, \hat{c}_2)$.} 
\label{fig9}
\end{figure}


\subsubsection*{Recovery of $H_2$ and $c_2$}
This statement also generates non convex functional, but now it has
excellent resolution with respect to $H_2$ 
(see Figures \ref{fig10} and \ref{fig11}).

\begin{figure}[ht]
\begin{center}
 \includegraphics[width=0.7\textwidth]{fig10} % surfW05005.jpg
\end{center}
\caption{Minimization of the functional $F(H_2, \hat{c}_2)$.}
\label{fig10}
\end{figure}

\begin{figure}[ht]
\begin{center}
 \includegraphics[width=0.7\textwidth]{fig11} % contour05005.jpg
\end{center}
\caption{Level line functional $F(H_2, \hat{c}_2)$.} 
\label{fig11}
\end{figure}

\subsection*{Conclusions}
In this publication we have shown how to derive mathematical models
for composite  media using its microstructure. As a rule, there is
some set of models depending on given criteria
$\mu_0,\,\lambda_0,\dots $ of the physical process in
consideration. For a fixed set of criteria the corresponding model
describes some of the main features of the process.

In the paper the simplest inverse problem was dealt with - recovery
of elastic parameters of the layer by Nelder-Mead algorithm. In the
future we are planning to establish connection upscaling procedure
and scattered waves and apply on this base recent developments of
true-amplitude imaging on the base of Gaussian beams for both
reflected and scattered waves \cite{Protasov_Tcheverda,Protasov_Tcheverda_1}.

\begin{thebibliography}{00}

\bibitem{Biot1} M. Biot;
\emph{Theory of propagation of elastic waves in a
fluid-saturated porous solid. I. Low-frequency range}. Journal of the
Acoustical Society of America. \textbf{28}, 168--178 (1955).

\bibitem{Biot2} M. Biot;
\emph{Theory of propagation of elastic waves in a
fluid-saturated porous solid.  II. Higher frequency range}. Journal
of the Acoustical Society of America. \textbf{28}, 179--191 (1955).

\bibitem{Biot3} M. Biot;
\emph{Generalized theory of seismic propagation in porous
dissipative media}. J. Acoust. Soc. Am. \textbf{34}, 1256--1264
(1962).

\bibitem{Berr1} J. G. Berryman;
\emph{Seismic wave attenuation in fluid-saturated porous media},  
J. Pure Appl. Geophys. (PAGEOPH) 128, 423-432 (1988).

\bibitem{Berr2} J. G. Berryman;
\emph{Effective medium theories for multicomponent
poroelastic composites}, ASCE Journal of Engineering Mechanics 132
(5), 519--531 (2006).

\bibitem{B-K} R. Burridge,  J. B. Keller;
\emph{Poroelasticity equations derived from microstructure}, 
 Journal of Acoustic Society of America  70,  No. 4  (1981),  1140 -- 1146.

\bibitem{S-P} E. Sanchez-Palencia;
\emph{Non-Homogeneous Media and Vibration Theory,
Lecture Notes in Physics},  Vol.  129, Springer-Verlag,  1980.

\bibitem{AM1} A. Meirmanov;
\emph{Mathematical models for poroelastic flows}, Atlantis Press, Paris, 2013.

\bibitem{AM2} A. Meirmanov;
\emph{A description of seismic seismic wave propagation in
porous media via homogenization}. SIAM J. Math. Anal. \textbf{40},
N3, 1272--1289 (2008).

\bibitem{AnM} A. Mikeli\'{c};
\emph{On the Efective Boundary Conditions Between Diferent Flow Regions.
Proceedings of the 1. Conference on Applied Mathematics and 
Computation Dubrovnik}, Croatia, September 13–18 (1999), 21-–37.

\bibitem{B-L-P} A. Bensoussan, J. L. Lions, G. Papanicolau;
 Asymptotic Analysis for Periodic Structures,
North-Holland, Amsterdam, 1978.

\bibitem{B-P} N. Bakhvalov, G. Panasenko;
\emph{Homogenization: averaging processes in periodic media},
 Math. Appl., Vol. 36, Kluwer Academic Publishers, Dordrecht, 1990.

\bibitem{GN1} G. Nguetseng;
\emph{A general convergence result for a
functional related to the theory of homogenization. }
 SIAM J. Math. Anal.  20 (1989), 608--623.

\bibitem{GN2} G. Nguetseng;
\emph{Asymptotic analysis for a stiff variational problem arising in mechanics}, 
 SIAM J. Math. Anal. 21 (1990), 1394--1414.

\bibitem{A-F} A. Fichtner;
\emph{Full Seismic Waveform Modelling and Inversion,
Springer- Verlag  Berlin Heidelberg}. (2011) pp. 113--140.

\bibitem{Protasov_Reshetova_Tcheverda_2}
M. I. Protasov, G. V. Reshetova, V. A. Tcheverda;
\emph{Fracture detection by Gaussian beam imaging of
seismic data and image spectrum analysis}, Geophysical Prospecting,
in-press (2015).

\bibitem{M-F} J. H. Mathews, K. D. Fink;
\emph{Numerical methods using Matlab}. 4th edition, 2004, Prentice-Hall Inc., 
chapter 8, pp.430--436.

\bibitem{Protasov_Tcheverda} M. I. Protasov, V. A. Tcheverda;
\emph{True amplitude  imaging by inverse generalized
Radon transform based on Gaussian beam decomposition of the acoustic
Green's function}, Geophysical Prospecting, \textbf{59}(2), 197--209 (2011).

\bibitem{Protasov_Tcheverda_1} M. I. Protasov, V. A. Tcheverda;
\emph{True-amplitude elastic  Gaussian beam imaging of multi-component walk-away VSP
data}, Geophysical Prospecting, \textbf{60}(6), 1030--1042 (2012).

\end{thebibliography}

\section*{Addendum posted on November 28, 2016}

The editor from  Zentralblatt informed us that a big portion of this article
coincides with the article 

``Seismic in composite media: elastic and poroelastic components''
by Anvarbek Meirmanov; Saltanbek Talapedenovich Mukhambetzhanov and
 Marat Nurtas (Sib. Elekron. Mat. Izv. 13, (2016) 75-88)  (Zbl 06607056).

The Electron. J. Differential Equations requested an explanation from the authors.
They replied that two co-authors submitted the manuscript to two 
different journals, and each eventually published it without consulting 
the other. They write,

\begin{quote}
It's my fault that I did not control the process. 
There is not any other explanation. 
Now I do not know what I should do. Maybe the best way here is to 
remove the paper from the site, if it is possible. 
I apologize once again,\\
yours sincerely,\\
Anvarbek Meirmanov.
\end{quote}

Since the article is already published, the EJDE editor has 
posted this explanation. We recommend that co-authors inform each other 
about their submissions. 
End of addendum.

\end{document}