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

\AtBeginDocument{{\noindent\small
2014 Madrid Conference on Applied Mathematics in honor of Alfonso Casal,\\
\emph{Electronic Journal of Differential Equations},
Conference 22 (2015),  pp. 1--18.\newline
ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu
\newline ftp ejde.math.txstate.edu}
\thanks{\copyright 2015 Texas State University.}
\vspace{9mm}}

\begin{document} \setcounter{page}{1}
\title[\hfilneg EJDE-2015/Conf/22 \hfil A viscoelastic-gravitational
layered earth model]
{Mathematical analysis of a viscoelastic-gravitational
layered earth model for magmatic intrusion in the dynamic case}

\author[A. Arjona, J. I. D\'iaz \hfil EJDE-2015/conf/22 \hfilneg]
{Alicia Arjona, Jes\'us Ildefonso D\'iaz}

\address{Alicia Arjona \newline
European Center for Geodynamics and Seismology,
Rue Josy Welter, 19, L-7256 Walferdange,
Gran-Duchy of Luxembourg}
\email{alicia.arj@gmail.com}

\address{Jes\'us Ildefonso D\'iaz \newline
Departamento de Matem\'atica Aplicada\\
Instituto de Matem\'atica Interdisciplinar\\
Universidad Complutense de Madrid,
Plaza de las Ciencias, 3, 28040 Madrid, Spain}
\email{jidiaz@ucm.es}

\thanks{Published November 20, 2015}
\subjclass[2010]{35K10, 35L10, 35Q86, 35Q74, 46E35, 86A60}
\keywords{Gravity changes; viscoelastic-gravitational earth model;
\hfill\break\indent weak solution; iterative algorithm; continuous dependence;
uniqueness of solutions}

\begin{abstract}
 Volcanic areas present a lower effective viscosity than usually in the
 Earth's crust. It makes necessary to consider inelastic properties in
 deformation modelling. As a continuation of work done previously by some of
 the authors, this work is concerned with the proof that the perturbed
 equations representing the viscoelastic-gravitational displacements
 resulting from body forces embedded in a layered Earth model leads to a
 well-posed problem even for any kind of domains, with the natural boundary
 and transmission conditions. A homogeneous or stratified viscoelastic
 half-space has often been used as a simple earth model to calculate the
 displacements and gravity changes. Here we give a constructive proof of
 the existence of weak solutions and we show the uniqueness and the continuous
 dependence with respect to the initial data of weak solutions of the dynamic coupled
 viscoelastic-gravitational field equations.
\end{abstract}

\maketitle
\numberwithin{equation}{section}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{remark}[theorem]{Remark}
\newtheorem{definition}[theorem]{Definition}
\allowdisplaybreaks


\section{Introduction} \label{intro}

The study of ground deformation in volcano has been an important issue
during last decades. There is a wide amount of literature on methodologies
for modelling elastic and viscoelastic response when a source is embedded in
media. The Mogi model \cite{Mogi1958} is the simplest analytical solution for a
point source of pressure in an elastic half-space to interpret ground
deformation. However, pure elastic models do not allow to reproduce gravity
changes in some events. Therefore, the computation of gravity changes and
deformations is advisable in order to do a correct interpretation. For an
example, see the following and references therein
\cite{Rundle1980,Rundle1981a,Rundle1982b,Rundle1983} and
\cite{Fernandez1994a,Fernandez1994b,Fernandez1999,Fernandez2001a,Fernandez2001b,
Fernandez2001c}.
Analytical and numerical solutions for modelling ground deformations and
gravity changes have been devised and used in literature
\cite{Arjona2011,Battaglia2004, Battaglia2008,Battaglia2009,
Camacho2011,Charco2014,Okada1958,Okubo1992,Yang1998}.
 These models consider different source geometries representing magma
such as spherical sources \cite{Mogi1958}, ellipsoidal point sources
\cite{Battaglia2004} and sources due to pressurisation of the magma
chamber \cite{Bonafede1998}. As we shall see the coupling with the
stationary equation for the potential gravity leads to many new difficulties
with respect to the pure viscoelastic models (see, e.g. \cite{Diaz2002}
and its references).

In volcanic areas the presence of inhomogeneous materials and high
temperature bodies reduce the effective viscosity of the Earth's crust.
Therefore, inelastic properties of the Earth's crust must be taken into
account \cite{Bonafede1998,Fernandez2001c}.
In this way Mogi's model was generalized to viscoelastic rheology
in \cite{Bonafede1998}. Rundle \cite{Rundle1982b} presented a stratified
viscoelastic half-space taking into account the interaction between the mass
of the intrusion and the ambient gravity field and the effect caused by the
change of pressure in the magmatic system. A theoretical and computational
methods for the calculation of viscoelastic-gravitational displacements
resulting from strike-slip faulting was described in \cite{Yu1996}.
The flow properties of the medium must also be considered. A Maxwell
viscoelastic fluid was used by Pollitz \cite{Pollitz1997}
 instead of Maxwell rheologies
\cite{Fernandez1996,Fernandez2001a,Fernandez2004,Folch2000,Rundle1982b}.
In this last case, the solution of the governing
equations can be obtained from elastic solution employing the correspondence
principle \cite{Fung1965}. A propagator matrix technique is used to obtain the
analytical solution of the elastic problem (see its description
in \cite{Fernandez1994a,Rundle1980}).

The objective of this work is to prove that the perturbed equations
representing the viscoelastic-gravitational displacements resulting from
body forces embedded in a layered Earth model leads to a well-posed problem
even for any kind of domains, with the natural boundary and transmission
conditions. The existence and uniqueness of weak solutions of the
elastic-gravitational problem was demonstrated in \cite{Arjona2008}
and the stabilization to solutions of the associated stationary
system was proved in \cite{Arjona2007}. We give here an additional
constructive proof of the existence of weak solutions and we show the
uniqueness and the continuous dependence with respect to the initial data of
weak solutions of the coupled viscoelastic-gravitational field equations.

\section{The problem} \label{prob}

We consider here an Earth model composed of several
viscoelastic-gravitational layers. We also consider the contribution of
source terms which represent magmatic intrusion, corresponding to body
forces acting on the medium. This is due to both volumetric change of wall
of the chamber and sudden emplacement of a mass into the medium as a result
of a new material injection into a magmatic chamber. The coupled model for
deformation and variation of gravity is given by the following system of
partial differential equations:
\begin{equation} \label{eP}
\begin{gathered}
\begin{aligned}
&\rho \mathbf{u( t,\mathbf{x})}_{tt}
 -\gamma \Delta \mathbf{u( t,\mathbf{x})}_t
 -\Delta \mathbf{u( t,\mathbf{x})}
 -\frac{1}{1-2\nu }\nabla ( \operatorname{div}
\mathbf{u( t,\mathbf{x})})
\\
&-\frac{\rho g}{\mu }\nabla ( \mathbf{u( t,\mathbf{x})}
\cdot \mathbf{e}_{z})
+ \frac{\rho g}{\mu }\mathbf{e}_{z}\operatorname{div}
 \mathbf{u( t,\mathbf{x})}+\frac{\rho }{\mu }\nabla
\phi( t,\mathbf{x})  \\
&=\mathbf{f}_{u}( t,\mathbf{x}),
\end{aligned}
\\
-\Delta \phi( t,\mathbf{x}) -4\pi
\rho G\operatorname{div}\mathbf{u( t,\mathbf{x})}
=f_{\phi }( t,\mathbf{x}) \quad \text{in  }(0,T)\times \Omega,
\end{gathered}
\end{equation}
where  $\mathbf{u}$ denotes the displacement, $\phi $ gravitational perturbed
potential, $\nu $ the Poisson's ratio, $\rho $ the unperturbed density of
the medium, $g$ the externally imposed gravitational acceleration, $\mu $ is
the rigidity, $\gamma \Delta \mathbf{u}_t$ is a term introduced by the
viscoelasticity of each layer, $G$ universal gravitational constant,
$\mathbf{e}_{z}$ is the unit vector pointing in the positive $z$-direction
(down into the medium) and $\mathbf{f}_{u}$ and $f_{\phi}$ the body forces.
Let us consider spatial domain $\Omega$ as union of $p$ layers ``overlay'', that
we will denote $\Omega _i$, $i=1,\dots,p$ . Each layer is given through
a common horizontal open set, $\omega \subset \mathbb{R}^2$, and so
\begin{equation}
\Omega_1:=\omega \times ( d_1,d_1+d_2),\quad
\Omega_2:=\omega \times ( d_1+d_2,d_1+d_2+d_{3}) , \dots;
\end{equation}
that is,
\begin{equation}
\Omega _i:=\omega \times \Big(
\sum_{j=1}^{i-1}d_{j},\sum_{j=1}^id_{j}\Big) \subset
\mathbb{R} ^{3},
\end{equation}
where $i=1,\dots,p-1$, and
\begin{equation}
\Omega _{p}:=\omega \times ( H,H+d_{r}) ,
\end{equation}
where $H:=\sum_{j=1}^{i-1}d_{j}$ and $d_{r}$ can be equal to $+\infty $.

Let $\mathbf{u}^i:[0,T]\times \Omega _i\to \mathbb{R}^{3}$
be the displacement vector in each layer where $T$ is an arbitrary
time, $\mathbf{u}^i=( u_{x}^i,u_{y}^i,u_{z}^i) $ which
depends on $\mathbf{x=}( x,y,z) $ and $t\in \left[ 0,T\right] $.
The system \eqref{eP} has been reached on each layer.

To the set of partial differential equations we will add the following
boundary conditions (see Figure \ref{fig1}). Regarding to displacement field we
prescribe on the side boundary, $\partial _{l}\Omega _i$, for $i=1,\dots,p $, that:
\begin{equation}
\mathbf{u}^i( t,\mathbf{x}) =\mathbf{0} \quad \mathbf{x}\in
\partial _{l}\Omega _i,t\in (0,T)
\label{cu1}
\end{equation}
on the upper boundary of the first layer $\partial _{+}\Omega _1$,
\begin{equation}
\frac{\partial \mathbf{u}^{1}( t,\mathbf{x}) }{\partial z}=
\mathbf{0} \quad \mathbf{x}\in \partial _{+}\Omega _1,t\in (0,T),
\label{cu2}
\end{equation}
and that on the bottom boundary $\partial _{-}\Omega _{p}$,
\begin{equation}
\mathbf{u}^p( t,\mathbf{x}) =\mathbf{0} \quad \mathbf{x}\in
\partial _{-}\Omega _{p},t\in (0,T).
\label{cu3}
\end{equation}

\begin{figure}[htb]
\begin{center}
\includegraphics[width=0.9\textwidth]{fig1}
\caption{Domain of the problem.}
\end{center}
\label{fig1}
\end{figure}

In general, we can assure only that the first derivatives of $\mathbf{u}$
are continuous on the boundaries of the layers, that is, on the boundary
between layers. We will require "transmission conditions" between both upper
and bottom boundaries of the layers excepting on the first and the last
layers. Therefore, the next conditions on
$\partial _{-}\Omega _i=\partial _{+}\Omega _{i+1}$, with
$i=1,\dots ,p-1$ are as follows:
\begin{equation}
\begin{gathered}
\mathbf{u}^i( t,\mathbf{x}) =\mathbf{u}^{i+1}( t,\mathbf{x}
) \quad \mathbf{x}\in \partial _{-}\Omega _i,\; t\in (0,T), \\
\frac{\partial \mathbf{u}^i( t,\mathbf{x}) }{\partial z}=
\frac{\partial \mathbf{u}^{i+1}( t,\mathbf{x}) }{\partial z} \quad
\mathbf{x}\in \partial _{-}\Omega _i,\; t\in (0,T).
\end{gathered}  \label{cu4}
\end{equation}

In relation to the gravitational perturbed potential we will assume that on side
boundary $\partial _{l}\Omega _i$ for $i=1,\dots,p$, it holds:
\begin{equation}
\phi ( t,\mathbf{x}) =0 \quad \mathbf{x}\in \partial _{l}\Omega
_i,\;t\in (0,T),
\label{cf1}
\end{equation}
on the upper boundary of the first layer $\partial _{+}\Omega _1$;
\begin{equation}
\phi ^{1}( t,\mathbf{x}) =\phi _0( t,\mathbf{x}) \quad
\mathbf{x}\in \partial _{+}\Omega _1,t\in (0,T),
\label{cf2}
\end{equation}
and on the bottom boundary, $\partial _{-}\Omega _{p}$;
\begin{equation}
\phi ^p( t,\mathbf{x}) =0 \quad \mathbf{x}\in \partial _{-}\Omega
_{p},t\in (0,T).
\label{cf3}
\end{equation}

Like before, we will require transmission conditions between upper and bottom
boundary of the next layers excepting on the first and the last
layers. So, we must have, on $\partial _{-}\Omega _i=\partial _{+}\Omega
_{i+1}$ with $i=1,\dots,p-1$, the following conditions:
\begin{equation}
\begin{gathered}
\phi ^i( t,\mathbf{x}) =\phi ^{i+1}( t,\mathbf{x}) \quad
\mathbf{x}\in \partial _{-}\Omega _i,t\in (0,T), \\
\frac{\partial \phi ^i( t,\mathbf{x}) }{\partial z}=\frac{
\partial \phi ^{i+1}( t,\mathbf{x}) }{\partial z} \quad \mathbf{x}\in
\partial _{-}\Omega _i,t\in (0,T).
\end{gathered}  \label{cf4}
\end{equation}

Finally we prescribe initial conditions for the displacements and the
velocities:
\begin{equation}
\begin{gathered}
\mathbf{u}( 0,\mathbf{x}) =\mathbf{u}_0( \mathbf{x})
\quad \text{in  }\Omega , \\
\mathbf{u}_t( 0,\mathbf{x}) =\mathbf{v}_0( \mathbf{x}
) \quad \text{in  }\Omega .
\end{gathered}  \label{ct}
\end{equation}

Since the multilayered structure of the domain introduce some possible
abrupt changes on the second derivatives of solutions the existence of
classical solutions of the problem looks artificial and we must introduce a
suitable weak formulation notion of the problem which mathematical analysis
is the main object of this paper.

\begin{remark} \label{rmk2.1}\rm
Here, and in what follows, $H^{1}(\Omega )$ denotes the Sobolev space
\begin{equation*}
H^{1}(\Omega )=\{\psi \in L^2(\Omega ):\frac{\partial \psi }{\partial
x_i}\in L^2(\Omega ),i=1,2,3\},
\end{equation*}
where $L^2$ is the space of all square integrable functions. Both spaces
have structure of Hilbert space (see, e.g. \cite{Brezis1984}, for more
details).
\end{remark}

\section{Weak formulation} \label{wf}

Following some ideas already introduced in the previous work by the authors
concerning the stationary problem \cite{Arjona2008} we define
the energy space $V=V_{u}\times V_{\phi }$ as cross product of the energy
spaces for the displacement and for the perturbed gravitational potential,
$V_{u}$ and $V_{\phi }$ respectively, where
\begin{align*}
V_{\phi }&:=\Big\{( \phi ^{1},\dots ,\phi ^p) \in
\prod_{i=1}^pH^{1}( \Omega _i) \text{: }\phi ^i=0
\text{ on  }\partial _{l}\Omega _i\forall \text{ }i=1\dots p,\phi
^{1}=0\text{ on }\partial _{+}\Omega _1\\
&\quad \text{and } \phi ^p=0\text{ on }\partial _{-}\Omega _{p},\phi ^i
=\phi ^{i+1},\; \frac{\partial \phi ^i( \mathbf{x}) }{\partial z}=
\frac{\partial \phi ^{i+1}( \mathbf{x}) }{\partial z}\\
&\quad \text{on } \partial _{-}\Omega _1\cup \partial _{+}\Omega _2\cup \partial
_{-}\Omega _2\cup \text{\dots }\cup \partial _{+}\Omega _{p},
\text{ for }i=1\dots p-1\Big\}
\end{align*}
and
\begin{align*}
V_{u}& :=\Big\{( \mathbf{u}^{1},\dots ,\mathbf{u}^p) \in
\prod_{i=1}^pH^{1}( \Omega _i) ^{3}:\mathbf{u}^i=0
\text{  on  }\partial _{l}\Omega _i\text{ for }i=1\dots p,
\mathbf{u}^{1}=0\text{  on  }\partial _{+}\Omega _1
\\
&\quad \text{and  }\mathbf{u}^p=0
\text{  on  }\partial _{-}\Omega _{p},\mathbf{u}^i=
\mathbf{u}^{i+1}\text{  and  }\frac{\partial \mathbf{u}
^i}{\partial z}=\frac{\partial \mathbf{u}^{i+1}}{\partial z}  \\
&\quad \text{on  }\partial _{-}\Omega _1\cup \partial _{+}\Omega _2\cup
\partial _{-}\Omega _2\cup \dots \cup \partial _{+}\Omega _{p},
\text{ for } i=1\dots p-1\Big\}.
\end{align*}

Regarding boundary data, $\phi _0$ shall be extended to the interior of
the domain $\Omega _1$ such that
\begin{equation}
\begin{gathered}
\widehat{\phi }_0\in L^{q}(0,T:H^{1}(\Omega _1)),\quad
\widehat{\phi } _0( t,\mathbf{x}) =\phi _0( t,\mathbf{x})
\text{ in }(0,T)\times \partial _{+}\Omega _1,\\
\widehat{\phi } _0( \mathbf{x}) =0
\text{ in } (0,T)\times (\partial _{-}\Omega _1\cup \partial _{l}\Omega _1),
\text{ for some }2\leq q\leq +\infty.
\end{gathered}  \label{fit}
\end{equation}
We will suppose, at least, that
\begin{equation}
\phi _0\in L^{q}(0,T:\prod_{i=1}^pH^{1}( \Omega _i) )
\text{ and satisfies \eqref{fit},}  \label{rdt1}
\end{equation}
and under the following regularity on the data:
\begin{gather}
\mathbf{f}_{u}\in L^2(0,T:\prod_{i=1}^pH^{-1}( \Omega _i) ^{3}),  \label{rdt2}\\
f_{\phi }\in L^{q}(0,T:\prod_{i=1}^pH^{-1}( \Omega
_i) )\text{, }  \label{rdt3}
\end{gather}
for some $2\leq q\leq +\infty $ and
\begin{equation*}
\mathbf{u}_0,\mathbf{v}_0\in V_{u}.
\end{equation*}

To introduce the weak solution definition we will follow similar arguments
already introduced in the previous work by the authors for the stationary
case. We start by assuming that $( \mathbf{u,}\phi ) $ is a
classical solution of  system \eqref{eP}. Let $( \mathbf{w,}\theta )
\in C^2([0,T]:V_{u}\times V_{\phi })$ be test functions. We multiply the
first equation of  \eqref{eP} by $\mathbf{w}^i( t,\mathbf{x}
) $, and the second equation by $\theta ^i( t,\mathbf{x})$.
Integrating by parts and applying Green's formula, we arrive, in a
natural way, to the definition of the weak solution of the problem.

\begin{definition} \label{def3.1} \rm
We assume the above regularity on the functions $\mathbf{f}_{u}$,
$f_{\phi }$, $\phi _0$, $\mathbf{u}_0$ and $\mathbf{v}_0$. We say that
$\{\mathbf{u},\phi\} $ is a weak solution of the problem \eqref{eP} with the
above mentioned boundary conditions   if
$( \mathbf{u,}\phi -\phi _0) \in L^2(0,T:V)$,
$\mathbf{u}_{tt}\in L^2(0,T:V_{u}')$ and for any test function
$( \mathbf{w}, \theta ) \in L^2(0,T:V)$,
$\mathbf{v\in }$ $H^{1}(0,T:V_{u}')$ the following equalities hold:
\begin{equation}
\begin{aligned}
&\int_0^{T}\sum_{i=1}^p\Big[ \langle \rho ^i
\mathbf{u}_{tt}^i,\mathbf{w}^i\rangle +\int_{\Omega _i}
\Big\{ \frac{1}{1-2\nu ^i}\operatorname{div}\mathbf{u}^i
\operatorname{div}\mathbf{w}^i
\\
&-\frac{\rho ^ig}{\mu ^i}\nabla ( \mathbf{u}^i\cdot \mathbf{
e}_{z}) \cdot \mathbf{w}^i
 +\frac{\rho ^ig}{\mu ^i}\mathbf{e}
_{z}\operatorname{div}\mathbf{u}^i \mathbf{w}^i+\nabla \mathbf{u}
^i:\nabla \mathbf{w}^i
+\gamma ^i\nabla \mathbf{u}
_t^i:\nabla \mathbf{w}^i\Big\} d\mathbf{x}\Big] dt
\\
&=\int_0^{T}\sum_{i=1}^p\Big[ \frac{\rho ^i}{\mu
^i}\int_{\Omega _i}-\nabla \phi ^i\cdot \mathbf{w}^i \, d
\mathbf{x}+\langle \mathbf{f}_{u}^i(t,\cdot),
\mathbf{w}^i(t,\cdot )\rangle _{V_{u}'\times V_{u}}\Big] dt,
\end{aligned}
\end{equation}
and a.e. $t\in (0,T)$,
\begin{equation}
\begin{aligned}
&\sum_{i=1}^p\int_{\Omega _i}\nabla \phi ^i( t,
\cdot ) \cdot \nabla \theta ^i( t,
\cdot ) d\mathbf{x}\\
&= \sum_{i=1}^p\Big[ 4\pi \rho ^iG
\int_{\Omega _i}\operatorname{div}\mathbf{u} ^i( t,\cdot ) \theta ^i( t,
\cdot ) d\mathbf{x}
+ \langle f_{\phi }^i( t,\cdot) ,\theta ^i( t,\cdot)
\rangle _{V_{\phi }'\times V_{\phi }}\Big].
\end{aligned}
\end{equation}
\end{definition}

We shall prove that problem \eqref{eP} is well-posed in the Hadamard sense.

\begin{theorem}\label{thm1}
(i) Assumed the regularity on the data $\mathbf{f}_{u}$, $f_{\phi }$,
$\phi _0$, $\mathbf{u}_0$ and $\mathbf{v}_0$ then there exists a unique weak
solution $\{ \mathbf{u},\phi\} $ of the problem \eqref{eP}.
Moreover $\mathbf{u}_t\in L^{\infty }(0,T:L^2(\Omega ))$,
 $\mathbf{u\in }L^{\infty }(0,T:V_{u})$, $\phi \in L^{q}(0,T:V_{\phi })$
and there exists a positive constant $C$ (depending on $T$, $\Omega _i$
and the constants $\rho ^i$, $\mu ^i$, $\nu ^i,\gamma ^i$ and $G$)
such that the following continuous dependence estimate holds
\begin{equation}
\begin{aligned}
&\sup_{t\in [ 0,T]}\sum_{i=1}^p\Big[
\int_{\Omega _i}| \mathbf{u}_t^i| ^2d\mathbf{x}\Big]
 +\int_0^{T}\int_{\Omega _i}| \nabla \mathbf{u}_t^i| ^2 \,d\mathbf{x}\,dt
+\sup_{t\in [0,T]}\sum_{i=1}^p\int_{\Omega _i}| \nabla \mathbf{u}
^i| ^2d\mathbf{x}  \\
&+\sup_{t\in [ 0,T]}\sum_{i=1}^p\int_{\Omega
_i}(\operatorname{div}\mathbf{u}^i)^2\,d\mathbf{x}
+\int_0^{T}\sum_{i=1}^p \int_{\Omega _i}| \nabla \phi ^i(t,\mathbf{x})|
^2\,d\mathbf{x} \\
&\leq C\Big[\int_0^{T}\sum_{i=1}^p\| \mathbf{f}_{u}^i(t,\cdot)\|
_{H^{-1}}^2dt+\int_0^{T}\sum_{i=1}^p\| \mathbf{f}_{\phi
}^i(t,\cdot)\| _{H^{-1}}^2dt \\
&\quad +\int_{\partial _{+}\Omega _1}| \phi _0(
\mathbf{s}) \mathbf{v}_0( \mathbf{s}) \cdot \mathbf{n}| ds
 +\sum_{i=1}^p\int_{\Omega _i}| \mathbf{v} _0^i( \mathbf{x}) | ^2d\mathbf{x}
\\
&\quad +\sum_{i=1}^p\int_{\Omega _i}| \nabla
\mathbf{v}_0^i( \mathbf{x}) | ^2d\mathbf{x}
+\sum_{i=1}^p\int_{\Omega _i}| \nabla \mathbf{u}
_0^i( \mathbf{x}) | ^2d\mathbf{x}
\\
&\quad +\sum_{i=1}^p\int_{\Omega _i}\operatorname{div}\mathbf{u}
_0^i( \mathbf{x}) ^2d\mathbf{x}+\int_0^{T}\int_{\partial
_{+}\Omega _1}\big| \phi ^0( t,\mathbf{s}) \frac{
\partial }{\partial \mathbf{n}}\phi ^0( t,\mathbf{s})
\big| ds\Big].
\end{aligned}
\label{cont depende estimate}
\end{equation}

(ii) If in addition we assume that
\begin{equation}
\begin{gathered}
\widehat{\phi }_0\in H^{1}(0,T:H^{1}(\Omega _1)),\quad
\phi _0\in H^{1}(0,T:\prod_{i=1}^pH^{1}( \Omega _i) ), \\
f_{\phi }\in H^{1}(0,T:\prod_{i=1}^pH^{-1}( \Omega_i) ),
\end{gathered}
\end{equation}
then $\ \phi \in H^{1}(0,T:V_{\phi })$ and we have the additional
continuous dependence estimate
\begin{align*}
&\sup_{t\in [ 0,T]}\sum_{i=1}^p\Big[
\int_{\Omega _i}| \mathbf{u}_t^i| ^2d\mathbf{x}\Big]
+\int_0^{T}\int_{\Omega _i}| \nabla \mathbf{u}
_t^i| ^2\,d\mathbf{x}\,dt
+\sup_{t\in [0,T]}\sum_{i=1}^p\int_{\Omega _i}| \nabla \mathbf{u}
^i| ^2d\mathbf{x} \\
&\quad +\sup_{t\in [ 0,T]}\sum_{i=1}^p\int_{\Omega
_i}(\operatorname{div}\mathbf{u}^i)^2\,d\mathbf{x}
+\sup_{t\in [ 0,T]}\sum_{i=1}^p\int_{\Omega _i}| \nabla
\phi ^i(t,\mathbf{x})| ^2\,d\mathbf{x}
 \\
&\leq C\Big[\int_0^{T}\sum_{i=1}^p\| \mathbf{f}_{u}^i(t,\cdot)\|
_{H^{-1}}^2dt
+\int_0^{T}\sum_{i=1}^p\| \mathbf{f}_{\phi}^i(t,\cdot)\|_{H^{-1}}^2dt\\
&\quad +\int_0^{T}\sum_{i=1}^p\| \frac{\partial }{\partial t}
(\mathbf{f}_{\phi }^i)(t,\cdot)\| _{H^{-1}}^2dt
 +\int_{\partial _{+}\Omega _1}| \phi _0(\mathbf{s})
 \mathbf{v}_0( \mathbf{s}) \cdot \mathbf{n} | ds
 +\sum_{i=1}^p\int_{\Omega _i}| \mathbf{v}_0^i( \mathbf{x}) | ^2d\mathbf{x}
\\
&\quad +\sum_{i=1}^p\int_{\Omega _i}| \nabla
\mathbf{v}_0^i( \mathbf{x}) | ^2d\mathbf{x}
+\sum_{i=1}^p\int_{\Omega _i}| \nabla \mathbf{u}
_0^i( \mathbf{x}) | ^2d\mathbf{x}
+\sum_{i=1}^p\int_{\Omega _i}\operatorname{div}\mathbf{u}_0^i( \mathbf{x
}) ^2d\mathbf{x}
\\
&\quad +\int_0^{T}\int_{\partial _{+}\Omega _1}| \phi
^0( t,\mathbf{s}) \frac{\partial }{\partial \mathbf{n}}\phi
^0( t,\mathbf{s}) | ds+\int_0^{T}\int_{\partial
_{+}\Omega _1} \big| \phi ^0( t,\mathbf{s}) \frac{
\partial ^2}{\partial t\partial \mathbf{n}}\phi ^0( t,\mathbf{s}
) \big| ds\Big],
\end{align*}
for a suitable positive constant $C$ (depending on $T$, $\Omega _i$ and
the constants $\rho ^i$, $\mu ^i$, $\nu ^i,\gamma ^i$ and $G$).
\end{theorem}

The existence of weak solutions will be proved by means of an iterative
method (which can be very useful to justify the convergence of some
numerical algorithms) without requiring any additional time regularity to
function $\phi $. This allows a great generality on the data. In the second
part we shall prove a stronger regularity on the weak solution by a
``cancellation method'' related to the time differentiability of function
$\phi $. We show that this holds under some slightly stronger regularity on
the data.

\section{Proof of Theorem \ref{thm1} part (i)}
\label{PT1}

We shall prove the existence of a weak solution by splitting it in several
steps. We first consider two different uncoupled problems: the first one
when displacements are known and the second one in which the gradient of the
gravitational potential is given.

\subsection{Uncoupled problem for the potential} ($\mathbf{u}$ is assumed to be
known)

We assume that $\mathbf{u}$ is known, with
\begin{equation}
\mathbf{u}^i\in H^{1}(0,T:H^{1}( \Omega _i) ).  \label{u dado}
\end{equation}
Let us consider the following problem, which we denotes as
 $P_1[\phi_0^{1},\mathbf{u}^i,\mathbf{f}_{\phi }^i]$, over the energy space
 $L^2(0,T:V_{\phi })$:
\begin{equation}
\begin{gathered}
-\Delta \phi ^i=4\pi \rho ^iG\operatorname{div}\mathbf{u}^i( t,\mathbf{x})
+f_{\phi }^i( t,\mathbf{x})  \quad \text{in }(0,T)\times \Omega
_i, \\
\phi ^i=0 \quad \text{on  }(0,T)\times \partial _{l}\Omega
_i\;\forall i=1,\dots ,p, \\
\phi ^i=\phi ^{i+1},\quad
\frac{\partial \phi ^i}{\partial z}=\frac{\partial \phi ^{i+1}}{\partial z}
\quad \text{on }(0,T)\times \partial _{-}\Omega _i=(0,T)\times \partial
_{+}\Omega _{i+1},\\
\forall i=1,\dots ,p-1,
\\
\phi ^{1}=\phi _0^{1} \quad \text{on }(0,T)\times \partial _{+}\Omega _1, \\
\phi ^p=0 \quad \text{on }(0,T)\times \partial _{-}\Omega _{p}.
\end{gathered}   \label{put}
\end{equation}

\begin{definition} \label{def4.1} \rm
Assumed the above regularity, a function $\phi $ is a weak solution of the
problem $P_1[\phi _0^{1},\mathbf{u}^i,\mathbf{f}_{\phi }^i]$ if
$\phi ^{\ast }:=\phi -\phi _0\in V_{\phi }$ and for any test function
$\theta \in V_{\phi }$, and a.e. $t\in (0,T)$ we have
\begin{equation}
\sum_{i=1}^p\int_{\Omega _i}\nabla \phi ^{\ast
i}\cdot \nabla \theta ^i\mathbf{\,}d\mathbf{x\,}
 = \sum_{i=1}^p4\pi \rho ^iG\int_{\Omega _i}(\operatorname{div}\mathbf{u}
^i)\,\theta ^i\,d\mathbf{x}+\langle f_{\phi }\mathbf{,}\,\theta
\rangle _{V_{\phi }'\times V_{\phi }}.
\end{equation}
\end{definition}

The following result was shown in the previous paper by the authors.

\begin{theorem}[\cite{Arjona2008}]
 Assuming the above regularity on the data $
\mathbf{u}^i$, $f_{\phi }$ and $\phi _0$, there exists a weak
solution, $\phi $, of the problem
\eqref{put}.
Moreover, if we denote the Poincar\'e's constant
on $\Omega _i$ by $C(\Omega _i)$ we have the estimate
\begin{align*}
\sum_{i=1}^p\int_{\Omega _i}| \nabla \phi
^i(t,\mathbf{x})| ^2\,d\mathbf{x}
& \leq  \sum_{i=1}^pC(\Omega _i)4\pi \rho ^iG\int_{\Omega
_i}| \nabla \mathbf{u}^i(t,\mathbf{x})| ^2d\mathbf{x}\\
&\quad +\sum_{i=1}^p\| \mathbf{f}_{\phi }^i(t,\cdot)\|_{H^{-1}}^2
+2\int_{\partial _{+}\Omega _1}\phi ^0( t,\mathbf{s})
\frac{\partial }{\partial \mathbf{n}}\phi ^0( t,\mathbf{s}) ds.
\end{align*}
\end{theorem}

\begin{remark} \label{rmk4.3} \rm
Note that in fact, since we assume that
$\mathbf{u}^i\in H^{1}(0,T:H^{1}( \Omega _i) )$ it is easy to see that
$f_{\phi } \in H^{1}(0,T:\prod_{i=1}^pH^{-1}( \Omega _i))$
implies a more regularity with respect to the variable
$t:\phi ^i\in H^{1}(0,T:H^{1}( \Omega _i) )$.
\end{remark}

\subsection{Uncoupled problem for the displacements}
($\phi $ is assumed to be known)

Now we assume given $\phi $ is given and
\begin{equation}
\phi ^i\in L^2(0,T:H^{1}( \Omega _i) ).  \label{dada fi}
\end{equation}
Let us consider the following problem $P_2[\phi ^i,\mathbf{f}_{u}^i]$
in $L^2(0,T:V_{u})$:
\begin{equation}
\begin{gathered}
\begin{aligned}
&\rho ^i\mathbf{u}_{tt}^i-\gamma ^i\Delta \mathbf{u}_t^i-\Delta
\mathbf{u}^i-\frac{1}{1-2\nu ^i}\nabla ( \operatorname{div}\mathbf{u}^i)  \\
&-\frac{\rho ^ig}{\mu ^i}\nabla ( \mathbf{u}^i\cdot \mathbf{e}
_{z}) +\frac{\rho ^ig}{\mu ^i}\mathbf{e}_{z}\operatorname{div}\mathbf{u}^i\\
&=- \frac{\rho _i}{\mu _i}\nabla \phi ^i+\mathbf{f}_{u}^i
\quad  \text{in }(0,T)\times \Omega _i, &
\end{aligned}
\\
\mathbf{u}^i( 0,\mathbf{x}) =\mathbf{u}_0^i( \mathbf{x}
) \quad \text{in }\Omega _i,   \\
\mathbf{u}_t^i( 0,\mathbf{x}) =\mathbf{v}_0^i(
\mathbf{x}) \quad \text{in }\Omega _i,  \\
\mathbf{u}^i=0\quad \text{ on } \partial _{l}\Omega_i,\quad
i=1,\dots ,p,   \\
\mathbf{u}^i=\mathbf{u}^{i+1},\quad
\frac{\partial \mathbf{u}^i}{ \partial z}
 =\frac{\partial \mathbf{u}^{i+1}}{\partial z}
\quad \text{on }\partial _{-}\Omega _i=\partial _{+}\Omega
_{i+1}, \;i=1,\dots ,p-1,
\\
\mathbf{u}^{1}=0\quad \text{on }\partial _{+}\Omega _1,   \\
\mathbf{u}^p=0\quad \text{on }\partial _{-}\Omega _{p}\,.
\end{gathered}   \label{pft}
\end{equation}

\begin{definition} \label{def4.4} \rm
We assume the above mentioned regularity on the data $f_{\phi }$ and
$\phi _0$ and initial data. The function $\mathbf{u}$ is a weak solution of
problem \eqref{pft}  if $\mathbf{u}\in
H^{1}(0,T:V_{u})$, $\mathbf{u}_{tt}\in L^2(0,T:V_{u}')$
and for any test function $\mathbf{w}$ such that
$\mathbf{w}\in H^{1}(0,T:V_{u})$ and $\mathbf{w}\in $ $H^2(0,T:V_{u}')$ we have
\begin{equation}
\begin{aligned}
&\sum_{i=1}^p\Big[ \int_0^{T}(\rho
^i\langle \mathbf{u}_{tt}^i,\mathbf{w}^i\rangle dt+\frac{1
}{1-2\nu ^i}\int_0^{T}\mathbf{\int_{\Omega _i}}\operatorname{div}\mathbf{u}^i\operatorname{div}
\mathbf{w}^i\,d\mathbf{x}\,dt   \\
& -\frac{\rho ^ig}{\mu ^i}\int_0^{T}\mathbf{
\int_{\Omega _i}}\nabla ( \mathbf{u}^i\cdot \mathbf{e}
_{z}) \cdot \mathbf{w}^i\,d\mathbf{x}\,dt+\frac{\rho ^ig}{\mu ^i}
\int_0^{T}\mathbf{\int_{\Omega _i}}e_{z}\operatorname{div}\mathbf{u}^i
\mathbf{w}^i\,d\mathbf{x}\,dt \\
& +\int_0^{T}\int_{\Omega _i}
\nabla \mathbf{u}^i:\nabla \mathbf{w}^i\,d\mathbf{x}\,dt
+ \int_0^{T}\int_{\Omega _i}\gamma ^i\nabla \mathbf{u}_t^i:\nabla
\mathbf{w}^id\mathbf{x}\Big]  \\
&=\sum_{i=1}^p\Big[ -\frac{\rho ^i}{\mu ^i}
\int_0^{T}\int_{\Omega _i}\nabla \phi ^i\cdot \mathbf{w}
^i\mathbf{\,}d\mathbf{x}+\int_0^{T}\langle \mathbf{f}_{u}^i(
t,\mathbf{x}) ,\mathbf{w}^i( t,\mathbf{x}) \rangle
_{V_{u}'\times V_{u}}\Big] .
\end{aligned}
\end{equation}
\end{definition}

The following result is a non difficult adaptation of a more general result
presented in \cite[Theorem 6.1, Chapter 3]{Duvaut1972}.

\begin{theorem} \label{thm4.5}
Assumed the regularity on $f_{\phi }$, the initial data and \eqref{dada fi}
there exists a weak solution, $\mathbf{u}$, of the problem \eqref{pft}.
\end{theorem}

\begin{proof}
We define the two bilinear forms
$a_{u}:V_{u}\times V_{u}\to \mathbb{R}$,
$a_{u_t}^{\ast }( \mathbf{u}_t\mathbf{,w}):V_{u}\times
V_{u}\to \mathbb{R}$ and the linear form $L_{u}:V_{u}
\to \mathbb{R}$ as follows:
\begin{gather}
\begin{aligned}
a_{u}( \mathbf{u,w})
&:= \sum_{i=1}^p\Big[ \int_{\Omega _i}
\Big\{ \frac{1}{1-2\nu ^i}\operatorname{div}\mathbf{u}^i\operatorname{div}
\mathbf{w}^i
-\frac{\rho ^ig}{\mu ^i}\nabla ( \mathbf{u}^i\cdot
\mathbf{e}_{z}) \cdot \mathbf{w}^i\\
&\quad +\frac{\rho ^ig}{\mu ^i}\mathbf{e}
_{z}\operatorname{div}\mathbf{u}^i\,\mathbf{w}^i\,+\nabla \mathbf{u}
^i:\nabla \mathbf{w}^i
 +\gamma ^i\nabla \mathbf{u}
_t^i:\nabla \mathbf{w}^i\,d\mathbf{x}\Big\} \Big] ,
\end{aligned} \label{a_u}
\\
a_{u_t}^{\ast }( \mathbf{u}_t\mathbf{,w})
:=\sum_{i=1}^p\mathbf{\int_{\Omega _i}}\gamma ^i\nabla \mathbf{u
}_t^i:\nabla \mathbf{w}^i\,d\mathbf{x},
\\
\langle L_{_{u}}(t),\mathbf{w}\rangle
:= \sum_{i=1}^p\Big[ \frac{\rho ^i}{\mu ^i}
\int_{\Omega _i}-\nabla \phi ^i(t,\text{\textperiodcentered })\cdot
\mathbf{w}^i\mathbf{\,}d\mathbf{x}+\langle \mathbf{f}_{u}^i(t,\text{
\textbf{\textperiodcentered }}),\mathbf{w}^i(t,\text{\textbf{
\textperiodcentered }})\rangle _{V_{u}'\times V_{u}}\Big] .
\end{gather}
That the bilinear form $a_{u}(\cdot,\cdot)$ \text{\textperiodcentered }
 is continuous and coercive,
and that the lineal form $L_{u}( t) $ is continuous was proved in
\cite[Theorem 3]{Arjona2008}. On the other hand, the
same type of arguments shows that the form
$a_{u_t}^{\ast }( \mathbf{u}_t\mathbf{,w}) $ is also continuous on
$V_{u}\times V_{u}$. Then,
we can apply the same arguments of the proof of
\cite[Theorem 6.1, Chapter 3]{Duvaut1972},
to the special case of no constraint, to the
problem
\begin{equation*}
\sum_{i=1}^p\langle \rho ^i\mathbf{u}_{tt}^i,\mathbf{w}
^i(t,\cdot )\rangle _{V_{u}'\times V_{u}}+a_{u}( \mathbf{u,w}) +a_{u_t}^{\ast }(
\mathbf{u}_t\mathbf{,w})
=\langle L_{_{u}}(t),\mathbf{w} \rangle _{V_{u}'\times V_{u}}
\end{equation*}
\end{proof}
(by using the transmission conditions) and we obtain the result.

\subsection{Coupled system}

To proof the existence and uniqueness of solutions of the
coupled system an iterative method will be used.
Firstly, we shall construct two sequences
$\{ \mathbf{u}^{n}( t,\mathbf{x})\} $ and
$\{ \phi ^{n}( \mathbf{x})\} $ as
follows. We start with $\phi ^0( \mathbf{x}) $ vector which has
initial data, $\phi _0( \mathbf{x}) $, as a first component and
rest of components $0$. With this vector and \eqref{pft}
 problem the unique vector
$\mathbf{u}^{1}( t,\mathbf{x}) $ is obtained.
Taking this vector as known, we solve \eqref{put}
problem to get the solution
associated to $\phi ^{1}$. In this way we build the next sequences (which
allow us to introduce some notations which will be used in what follows):
\begin{equation*}
\phi ^0=\begin{pmatrix}
\phi _0 \\
0 \\
\cdot \\
0
\end{pmatrix}
 \overset{P_2[\phi _0^i,\mathbf{f}_{u}^i] }{\to }
\mathbf{u}^{1}
=\begin{pmatrix}
\mathbf{u}_1^{1} \\
\mathbf{u}_2^{1} \\
\cdot \\
\mathbf{u}_{p}^{1}
\end{pmatrix}
\overset{P_1[\phi _0^{1},\mathbf{u}^i,f_{\phi }^i]}{\to }\phi ^{1}
=\begin{pmatrix}
\phi _1^{1} \\
\phi _2^{1} \\
\cdot \\
\phi _{p}^{1}
\end{pmatrix}.
\end{equation*}
In general,
\begin{equation*}
\phi ^{n-1}=\begin{pmatrix}
\phi _1^{n-1} \\
\phi _2^{n-1} \\
\cdot \\
\phi _{p}^{n-1}
\end{pmatrix}
\overset{P_2[\phi ^i,\mathbf{f}_{u}^i] }{\to }
\mathbf{u}^{n}=\begin{pmatrix}
\mathbf{u}_1^{n} \\
\mathbf{u}_2^{1} \\
\cdot \\
\mathbf{u}_{p}^{n}
\end{pmatrix}
\overset{P_1[\phi _0^{1},\mathbf{u}^i,f_{\phi }^i] }{ \to }\phi ^{n}
=\begin{pmatrix}
\phi _1^{n} \\
\phi _2^{n} \\
\cdot \\
\phi _{p}^{n}
\end{pmatrix}.
\end{equation*}

We claim that it is possible to find some universal a priori
estimates (i.e., independent on $n$) allowing to pass to the limit.
Indeed, by multiplying the equation of $\mathbf{u}_i^{n}$ by
$\rho ^i\mathbf{u}_{i,t}^{n}( t,\mathbf{x}) $
(where $\mathbf{u}_{i,t}^{n}( t,\mathbf{x})
=\frac{\partial \mathbf{u}_i^{n}( t, \mathbf{x}) }{\partial t}$).
Since, in general, we know that
\begin{gather*}
\nabla \mathbf{u}^i\cdot \nabla \mathbf{u}_t^i
=\frac{1}{2}\frac{\partial }{\partial t}( | \nabla \mathbf{u}^i| ^2),
\\
\int_0^{T}\int_{\Omega _i}\Delta \mathbf{u}^i\cdot \mathbf{u}_t^i
=- \frac{1}{2}\int_0^{T}\frac{d}{dt}\int_{\Omega _i}| \nabla
\mathbf{u}^i| ^2.
\\
\int_0^{T}\int_{\Omega _i}\operatorname{div}\mathbf{u}^i
\operatorname{div}\mathbf{u}_t^i=
\frac{1}{2}\int_0^{T}\frac{d}{dt}\int_{\Omega _i}( \operatorname{div}\mathbf{u}
^i) ^2\,.
\end{gather*}
Then, by using Green's formula, denoting the Poincar\'{e}'s constant on
$\Omega _i$, by $C(\Omega _i)$, and applying Poincar\'{e}'s and Young's
inequalities ($ab\leq \varepsilon a^2+C_{\varepsilon }b^2$) we obtain
the estimate
\begin{equation}
\begin{aligned}
&\sup_{t\in [ 0,T]}\sum_{i=1}^p\rho ^i(\frac{1}{
2}-\frac{Tg\rho ^i}{\mu ^i})\int_{\Omega _i}| \mathbf{u}
_{i,t}^{n}| ^2d\mathbf{x}
+\sum_{i=1}^p(1-\varepsilon )\rho ^i\gamma
^i\int_0^{T}\int_{\Omega _i}| \nabla \mathbf{u}
_{i,t}^{n}| ^2\,d\mathbf{x}\,dt \\
&+\sup_{t\in [ 0,T]}\sum_{i=1}^p\frac{\rho ^i}{2
}(1-\frac{Tg\rho ^i}{\mu ^i})\int_{\Omega _i}| \nabla
\mathbf{u}_i^{n}| ^2d\mathbf{x}\\
&+\sup_{t\in [0,T]}\sum_{i=1}^p\frac{\rho ^i}{2}\big(\frac{1}{(1-2\nu ^i)}-
\frac{Tg\rho ^i}{\mu ^i}\big)
\int_{\Omega _i}(\operatorname{div}\mathbf{u}_i^{n})^2\,d\mathbf{x} \\
&\leq C_{\varepsilon }\sum_{i=1}^p\int_0^{T}\Big(\int_{\Omega _i}| \phi
_i^{n-1}(t,\mathbf{x})| ^2d\mathbf{x} +\frac{\rho ^i}{\gamma
^i}\| \mathbf{f}_{\mathbf{u}}^i(t,\cdot)\| _{H^{-1}}^2\Big)dt
\\
&\quad +\sum_{i=1}^p\frac{\rho ^i}{2}\int_{\Omega
_i}| \mathbf{v}_0^i( \mathbf{x}) | ^2d
\mathbf{x}+\sum_{i=1}^p\gamma ^i\rho ^i\int_{\Omega
_i}| \nabla \mathbf{v}_0^i( \mathbf{x})| ^2d\mathbf{x} \\
&\quad +\sum_{i=1}^p\frac{\rho ^i}{2}\int_{\Omega
_i}| \nabla \mathbf{u}_0^i( \mathbf{x})
| ^2d\mathbf{x}+\sum_{i=1}^p\frac{\rho ^i}{2(1-2\nu
^i)}\int_{\Omega _i}\operatorname{div}\mathbf{u}_0^i( \mathbf{x}) ^2d
\mathbf{x} \\
&\quad +\int_0^{T}\int_{\partial _{+}\Omega _1}| \phi
^0( t,\mathbf{s}) \frac{\partial }{\partial \mathbf{n}}\phi
^0( t,\mathbf{s}) | ds.
\end{aligned}
\label{Estimate iteration u}
\end{equation}
In the above estimate we used the following inequalitites
\begin{gather*}
\begin{aligned}
&\sum_{i=1}^p\frac{\rho ^i}{(\mu ^i)^2\gamma ^i}
\int_0^{T}\int_{\Omega _i}\nabla \phi _i^{n-1}(t,\mathbf{x})\cdot
\mathbf{u}_{i,t}^{n}( t,\mathbf{x}) d\mathbf{x}\,dt\\
&= -\sum_{i=1}^p\frac{\rho ^i}{(\mu ^i)^2\gamma ^i}
\int_0^{T}\int_{\Omega _i}\phi _i^{n-1}(t,\mathbf{x})\operatorname{div}\mathbf{u}
_{i,t}^{n}( t,\mathbf{x}) d\mathbf{x}\,dt \\
&\leq \sum_{i=1}^p(C_{\varepsilon }\int_0^{T}\int_{\Omega
_i}| \phi _i^{n-1}(t,\mathbf{x})| ^2d\mathbf{x}\,dt
+ \epsilon \rho ^i\gamma ^i\int_0^{T}\int_{\Omega _i}|
\nabla \mathbf{u}_{i,t}^{n}| ^2\,d\mathbf{x}\,dt),
\end{aligned}
\\
\big| \int_0^{T}\int_{\Omega _i}\nabla ( \mathbf{u}
_i^{n}\cdot \mathbf{e}_{z}) \cdot \mathbf{u}_{i,t}^{n}\big|
\leq \frac{T}{2}\sup_{t\in [ 0,T]}\int_{\Omega _i}| \nabla
\mathbf{u}_i^{n}| ^2+\frac{T}{2}\sup_{t\in [
0,T]}\int_{\Omega _i}| \mathbf{u}_{i,t}^{n}| ^2,
\\
\big| \int_0^{T}\int_{\Omega _i}\mathbf{e}_{z}\operatorname{div}\mathbf{u}
_i^{n}\cdot \mathbf{u}_{i,t}^{n}\big|
\leq \frac{T}{2}\sup_{t\in [ 0,T]}\int_{\Omega _i}|
\operatorname{div}\mathbf{u}_i^{n}|^2
+\frac{T}{2}\sup_{t\in [ 0,T]}\int_{\Omega _i}|
\mathbf{u}_{i,t}^{n}| ^2.
\end{gather*}

 On the other hand, since $L^2(\Omega _i)\subset H^{-1}(\Omega_i)$,
there exists $K>0$ such that
\begin{equation*}
\| \operatorname{div}\mathbf{u}_i^{n-1}\| _{H^{-1}(\Omega _i)}\leq
K\| \operatorname{div}\mathbf{u}_i^{n-1}\| _{L^2(\Omega _i)}\mathbf{.
}
\end{equation*}
Then, from the coerciveness of the bilinear form associated to the elliptic
equation we know that
\begin{equation}
\begin{aligned}
 \int_{\Omega _i}| \nabla \phi _i^{n}( t,\mathbf{x})
| ^2\,d\mathbf{x}
&\leq  K\Big(\int_{\Omega _i}| \operatorname{div}
\mathbf{u}_i^{n-1}( t,\mathbf{x}) | ^2\,d\mathbf{x}
+\sum_{i=1}^p\| \mathbf{f}_{\phi }^i(t,\cdot)\|_{H^{-1}}^2 \\
&\quad +\int_{\partial _{+}\Omega _1}
\big| \phi ^0( t,\mathbf{s}) \frac{\partial }{\partial \mathbf{n}}
\phi ^0( t,\mathbf{s}) \big| ds\Big).
\end{aligned}\label{estimate gradient fi}
\end{equation}
Using twice Poincar\'{e}'s inequality we obtain
\begin{align*}
&\sum_{i=1}^p\int_0^{T}\int_{\Omega _i}|
\phi _i^{n}(t,\mathbf{x})| ^2d\mathbf{x}\,dt\\
&\leq \sum_{i=1}^p\widehat{K}(C(\Omega _i))\int_0^{T}\int_{\Omega
_i}| \nabla \phi _i^{n}( t,\mathbf{x}) | ^2\,d\mathbf{x}\,dt
 \\
&\leq \sum_{i=1}^p\widehat{K}\Big\{
\int_0^{T}\int_{\Omega _i}| \mathbf{u}_i^{n-1}( T,
\mathbf{x}) | ^2\,d\mathbf{x}\,dt
+ \int_0^{T}\sum_{i=1}^p\| \mathbf{f}_{\phi }^i(t,\cdot)\| _{H^{-1}}^2dt \\
&\quad  +\int_{\partial _{+}\Omega _1}\big| \phi
^0( t,\mathbf{s}) \frac{\partial }{\partial \mathbf{n}}\phi
^0( t,\mathbf{s}) \big| ds\Big\}  \\
&\leq \sum_{i=1}^p\overline{K}
\Big\{ T\sup_{t\in [ 0,T]}\int_{\Omega _i}| \nabla \mathbf{u}_i^{n-1}(
t,\mathbf{x}) | ^2\,d\mathbf{x}
+\int_0^{T}\sum _{i=1}^p\| \mathbf{f}_{\phi }^i(t,\cdot)\|
_{H^{-1}}^2dt  \\
&\quad  +\int_{\partial _{+}\Omega _1}\big| \phi
^0( t,\mathbf{s}) \frac{\partial }{\partial \mathbf{n}}\phi
^0( t,\mathbf{s}) \big| ds\Big\} ,
\end{align*}
where $\widehat{K}$ and $\overline{K}$ are  positive constants depending
increasingly on the Poincar\'{e}'s constants $C(\Omega _i)$. Thus,
combining both inequalities we obtain
\begin{align}
&\sup_{t\in [ 0,T]}\sum_{i=1}^p
\rho ^i\big(\frac{1}{2}-\frac{Tg\rho ^i}{\mu ^i}\big)
\int_{\Omega _i}| \mathbf{u}_t^i| ^2d\mathbf{x}  \nonumber \\
& +\sum_{i=1}^p(1-\varepsilon C(\Omega _i))\rho
^i\gamma ^i\int_0^{T}\int_{\Omega _i}| \nabla
\mathbf{u}_t^i| ^2\,d\mathbf{x}\,dt  \nonumber\\
& +\sup_{t\in [ 0,T]}\sum_{i=1}^p
\frac{\rho ^i}{2}\big(1-\frac{Tg\rho ^i}{\mu ^i}
\big)\int_{\Omega _i}| \nabla \mathbf{u}^i| ^2d\mathbf{x} \nonumber \\
& +\sup_{t\in [ 0,T]}\sum_{i=1}^p\frac{\rho ^i}{2}\big(
\frac{1}{(1-2\nu ^i)}-\frac{Tg\rho ^i}{\mu ^i}\big)
\int_{\Omega _i}(\operatorname{div}
\mathbf{u}^i)^2\,d\mathbf{x} +\sum_{i=1}^p\int_0^{T}\int_{
\Omega _i}| \phi _i^{n}(t,\mathbf{x})| ^2d\mathbf{x} dt \nonumber \\
&\leq \sum_{i=1}^p\overline{K}\Big\{T\sup_{t\in [
0,T]}\int_{\Omega _i}| \nabla \mathbf{u}_i^{n-1}( t,
\mathbf{x}) | ^2\,d\mathbf{x} +\sum_{i=1}^p
\int_0^{T}\sum_{i=1}^p\| \mathbf{f}_{\phi
}^i(t,\cdot)\| _{H^{-1}}^2dt  \nonumber\\
&\quad +\sum_{i=1}^p\int_0^{T}\| \mathbf{f}_{\mathbf{u}
}^i(t,\cdot)\| _{H^{-1}}^2dt++\sum_{i=1}^p\int_{\Omega
_i}| \mathbf{v}_0^i( \mathbf{x}) | ^2d
\mathbf{x} +\sum_{i=1}^p\int_{\Omega _i}| \nabla \mathbf{v
}_0^i( \mathbf{x}) | ^2d\mathbf{x}
\nonumber \\
&\quad +\sum_{i=1}^p\int_{\Omega _i}| \nabla \mathbf{u}
_0^i( \mathbf{x}) | ^2d\mathbf{x}
\sum_{i=1}^p\int_{\Omega _i}\operatorname{div}\mathbf{u}_0^i( \mathbf{x
}) ^2d\mathbf{x}\\
&\quad +\int_0^{T}\int_{\partial _{+}\Omega _1}\big| \phi ^0( t,\mathbf{s})
 \frac{\partial }{\partial \mathbf{n}}\phi ^0( t,\mathbf{s}) \big| ds\,dt\Big\}.
\label{almost final ineq}
\end{align}
Summarizing, if we assume $T>0$ and $\varepsilon $ small enough we conclude
that any uniform estimate on the term
\begin{equation*}
I_{n}=\sup_{t\in [ 0,T]}\sum_{i=1}^p\int_{\Omega
_i}| \nabla \mathbf{u}_i^{n}( t,\mathbf{x})| ^2\,d\mathbf{x}
\end{equation*}
allows us to make uniform all the above inequalities.
But we can understood \eqref{almost final ineq} in the form
\begin{equation*}
I_{n}\leq \delta I_{n-1}+A
\end{equation*}
with
\begin{equation*}
\delta :=\frac{\overline{K}T}{\min_i\frac{\rho ^i}{2}\big(1-\frac{Tg\rho
^i}{\mu ^i}\big)}
\end{equation*}
and $A$ given trough the external and initial data.
Thus, if $\delta \in (0,1)$, we obtain
\begin{equation*}
\lim_{n}I_{n}\leq A
\end{equation*}
and we obtain the uniform estimate $I_{n}\leq 1+A$ for any
 $n\in \mathbb{N}$
large enough (which implies uniform estimates in \eqref{Estimate iteration u}
 and \eqref{estimate gradient fi}). But, as before, the condition
$\delta \in (0,1)$ holds if we assume $T=T_0>0$ small enough and thus we obtain a set
of uniform a priori estimates on the sequences
$\{ \mathbf{u}^{n}\} $ and $\{ \phi ^{n}\} $ which show that they
converge weakly in $H^{1}(0,T_0:V_{u})$ and $L^2(0,T_0:V_{\phi })$,
respectively, to a vectorial function $(\mathbf{u,}\phi )$ which is a local
(in time) weak solution of the coupled system.

To prove the uniqueness of the local weak solution (i.e. when
$t\in [ 0,T_0]$ it suffices to show that if we assume as zero
all the data then the unique solution is the trivial solution $(\mathbf{u,}
\phi )=(\mathbf{0},0)$. This is an special conclusion of the obtained
continuous dependence estimate.

 Moreover, such local time $T_0$ does not depend on any norm of
the data (but only on the coefficients). In particular, if $T=T_0$ we obtain
the estimate
\begin{equation}
\begin{aligned}
&\sum_{i=1}^p\int_0^{T}\int_{\Omega _i}|
\nabla \mathbf{u}_t^i| ^2\,d\mathbf{x}\,dt
+\sum _{i=1}^p\int_0^{T}\int_{\Omega _i}| \phi _i^{n}(t,
\mathbf{x})| ^2d\mathbf{x}\,dt \\
&\leq  \sum_{i=1}^p\overline{K}\{\sum_{i=1}^p\int_0^{T}
\sum_{i=1}^p\| \mathbf{f}_{\phi }^i(t,\cdot)\|_{H^{-1}}^2dt
 + \sum_{i=1}^p\int_0^{T}\| \mathbf{f}_{\mathbf{u}
}^i(t,\cdot)\| _{H^{-1}}^2dt\\
&\quad +\sum_{i=1}^p\int_{\Omega
_i}| \mathbf{v}_0^i( \mathbf{x}) | ^2d \mathbf{x}
+\sum_{i=1}^p\int_{\Omega _i}| \nabla \mathbf{v
}_0^i( \mathbf{x}) | ^2d\mathbf{x}
 \\
&\quad + \sum_{i=1}^p\int_{\Omega _i}| \nabla \mathbf{u}
_0^i( \mathbf{x}) | ^2d\mathbf{x}
+\sum_{i=1}^p\int_{\Omega _i}\operatorname{div}\mathbf{u}_0^i( \mathbf{x
}) ^2d\mathbf{x}\\
&\quad +\int_0^{T}\int_{\partial _{+}\Omega
_1}\big| \phi ^0( t,\mathbf{s}) \frac{\partial }{
\partial \mathbf{n}}\phi ^0( t,\mathbf{s}) \big| \,ds\,dt\},
\end{aligned}
\label{universal estimate}
\end{equation}
with $\overline{K}>0$ independent of $T_0$.
Thus we can iterate the estimate on the intervals $[mT_0,(m+1)T_0]$ and
the estimate \eqref{universal estimate} remains valid.
So, no possible blow-up of the norms of
$(\mathbf{u,}\phi )$ on $(0,T)$ may arise if $T>0$ is any arbitrary number
and part i) of Theorem 1 is completely proved (note that the resultant
constant $C$ in estimate \eqref{cont depende estimate} can be taken as
$C=\overline{K}/[1-cT_0]^{m}$ if $T\in [ mT_0,(m+1)T_0]$ for some
natural number $m$ and a suitable constant $c$ such that $cT_0<1$.

\section{Proof of Theorem \ref{thm1} part (ii)}
\label{PT12}

 In contrast with the above set of estimates we can use now
different arguments since we already know that there exists a solution
of the iterative algorithm. We use now an idea which we can call as a
``cancellation argument'' (among the equations and boundary conditions) which
essentially consists in differentiating the equation of $\phi ^i$ with
respect to $t$. If we neglect, for a while, the external data then we obtain
\begin{equation*}
\frac{\partial }{\partial t}( -\Delta \phi ^i)
=4\pi \rho ^iG \frac{\partial }{\partial t}\operatorname{div}\mathbf{u}^i
=4\pi \rho ^iG\operatorname{div}\mathbf{u} _t^i
\end{equation*}
and since the right hand side was controlled in the previous estimate we obtain
that the left hand side is integrable in a suitable functional space. More
precisely, by multiplying last expression by
$\frac{\rho ^i}{\mu ^i}\phi^i$ and integrating over the space we obtain
\begin{equation}
\int_{\Omega _i}\frac{\rho ^i}{\mu ^i}\nabla \phi
^i\cdot \nabla \phi _t^i~d\mathbf{x}=\int_{\Omega _i}
\frac{4\pi ( \rho ^i) ^2}{\mu ^i}G\operatorname{div}\mathbf{u}
_t^i\,\phi ^i\,d\mathbf{x}.
\label{eq derivate in t}
\end{equation}
But the right hand side of \eqref{eq derivate in t} arises also when we
multiply the equation of $\mathbf{u}_i$ by
$\rho ^i\mathbf{u} _{i,t}( t,\mathbf{x}) $.
Finally, if we apply such a process but
now taking into account the contributions of the body forces $f_{\phi }^i$
and $\mathbf{f}_{u}^i$, the ones of the boundary data (when integrating by
parts, specially on $\partial _{+}\Omega _1$), and the one of the initial
data we obtain
\begin{align*}
&\sup_{t\in [ 0,T]}\sum_{i=1}^p\Big[
\int_{\Omega _i}\Big( 2\pi ( \rho ^i) ^2G-\frac{T4\pi
( \rho ^i) ^2Gg}{\mu ^i}\Big) | \mathbf{u}_t^i| ^2d\mathbf{x}\Big]  \\
&+\int_0^{T}\int_{\Omega _i}4\pi \rho ^iG\gamma
^i| \nabla \mathbf{u}_t^i| ^2\,d\mathbf{x}\,dt
 \\
&+\sup_{t\in [ 0,T]}\sum_{i=1}^p\Big[
\int_{\Omega _i}\Big( 2\pi \rho ^iG-\frac{T2\pi ( \rho
^i) ^2Gg}{\mu ^i}\Big) | \nabla \mathbf{u}
^i| ^2d\mathbf{x}  \\
&+\Big( \frac{2\pi \rho ^iG}{1-2\nu ^i}-\frac{T2\pi
( \rho ^i) ^2Gg}{\mu ^i}\Big)
\int_{\Omega _i}(\operatorname{div}\mathbf{u}^i)^2\,d\mathbf{x} \Big]
\\
&\leq \sup_{t\in [ 0,T]}\sum_{i=1}^p\int_{\Omega_i}
 \frac{\rho ^i}{2\mu ^i}| \nabla \phi ^i(t,\mathbf{x} )| ^2\,d\mathbf{x} \\
&\quad +\sum_{i=1}^p4\pi \rho ^iG\rho
^i\int_0^{T}\langle \mathbf{f}_{u}^i,\mathbf{u}
_t^i\rangle dt+\sum_{i=1}^p\frac{\rho ^i}{\mu ^i}
\int_0^{T}\langle \frac{\partial }{\partial t}(f_{\phi }^i),\phi
^i\rangle dt \\
&\quad +\frac{4\pi ( \rho ^{1}) ^2G}{\mu ^{1}}
\int_{\partial _{+}\Omega _1}| \phi _0( \mathbf{s})
\mathbf{v}_0( \mathbf{s}) \cdot \mathbf{n}|
ds+\sum_{i=1}^p2\pi ( \rho ^i) ^2G\int_{\Omega
_i}| \mathbf{v}_0^i( \mathbf{x}) | ^2d\mathbf{x} \\
&\quad +\sum_{i=1}^p4\pi \rho ^iG\gamma ^i\int_{\Omega
_i}| \nabla \mathbf{v}_0^i( \mathbf{x})
| ^2d\mathbf{x}+\sum_{i=1}^p2\pi \rho ^iG\int_{\Omega
_i}| \nabla \mathbf{u}_0^i( \mathbf{x})
| ^2d\mathbf{x} \\
&\quad +\sum_{i=1}^p\frac{2\pi \rho ^iG}{1-2\nu ^i}
\int_{\Omega _i}\operatorname{div}\mathbf{u}_0^i( \mathbf{x}) ^2d\mathbf{
x}+\frac{\rho ^{1}}{\mu ^{1}}\int_0^{T}\int_{\partial _{+}\Omega
_1}\big| \phi ^0( t,\mathbf{s}) \frac{\partial ^2}{
\partial t\partial \mathbf{n}}\phi ^0( t,\mathbf{s})
\big| ds,
\end{align*}
where we used that
\begin{equation*}
\int_0^{T}\int_{\Omega _i}\frac{\rho ^i}{\mu ^i}\nabla \phi
^i\cdot \nabla \phi _t^i=\frac{1}{2}\int_0^{T}\frac{d}{dt}
\int_{\Omega _i}\frac{\rho ^i}{\mu ^i}| \nabla \phi
^i| ^2.
\end{equation*}
Similarly, for $t\in [ 0,T]$, using the Poincar\'{e}'s constants
$C(\Omega _i)$
\begin{align*}
&\sum_{i=1}^p\int_{\Omega _i}| \nabla \phi
^i(t,\mathbf{x})| ^2\,d\mathbf{x} \\
& \leq \sum_{i=1}^pC(\Omega _i)4\pi \rho ^iG\int_{\Omega
_i}| \nabla \mathbf{u}^i(t,\mathbf{x})| ^2d\mathbf{
x+}\sum_{i=1}^p\| \mathbf{f}_{\phi }^i(t,\cdot)\| _{H^{-1}}^2 \\
&\quad +2\int_{\partial _{+}\Omega _1}\big| \phi
^0( t,\mathbf{s}) \frac{\partial }{\partial \mathbf{n}}\phi
^0( t,\mathbf{s}) \big| ds.
\end{align*}
Adding both inequalities and by taking a suitable constant $K$
(depending increasingly on the Poincar\'{e}'s constant $C(\Omega _i)$)
we obtain the result.

\section{Associated stationary system}
\label{OSS}

The above arguments can be also applied to prove the uniqueness of the
associated stationary problem
\begin{equation} \label{Pinf}
\begin{gathered}
\begin{aligned}
&\rho ^i\mathbf{u}_{tt}^i-\gamma ^i\Delta \mathbf{u}_t^i-\Delta
\mathbf{u}^i-\frac{1}{1-2\nu ^i}\nabla ( \operatorname{div}\mathbf{u}^i)
-\frac{\rho ^ig}{\mu ^i}\nabla ( \mathbf{u}^i\cdot \mathbf{e}
_{z}) +\frac{\rho ^ig}{\mu ^i}\mathbf{e}_{z}\operatorname{div}\mathbf{u}^i\\
&=- \frac{\rho _i}{\mu _i}\nabla \phi ^i+\mathbf{f}_{u}( \mathbf{x})
\quad \text{in  }\Omega _i,
\end{aligned}
\\
-\Delta \phi ^i=4\pi \rho ^iG\operatorname{div}\mathbf{u}^i( t,\mathbf{x})
+f_{\phi }( \mathbf{x})  \quad \text{in    }\Omega_i,
\end{gathered}
\end{equation}
under the same type of boundary conditions (but now with data independent of
$t$):
\begin{equation} \label{BCinf}
\begin{gathered}
\mathbf{u}^i=0,\quad \phi ^i=0, \quad \text{on  }
\partial_{l}\Omega _i\; \forall i=1,\dots ,p, \\
\mathbf{u}^i=\mathbf{u}^{i+1}, \quad
\frac{\partial \mathbf{u}^i}{\partial z}=\frac{\partial \mathbf{u}^{i+1}}{\partial z}
\quad \text{on }\partial _{-}\Omega _i=\partial _{+}\Omega _{i+1},\;
i=1,\dots ,p-1,
 \\
\phi ^i=\phi ^{i+1},\quad \frac{\partial \phi ^i}{\partial z}
=\frac{\partial \phi ^{i+1}}{\partial z}
\quad \text{on }\partial _{-}\Omega _i=\partial _{+}\Omega _{i+1},\;
i=1,\dots ,p-1,
 \\
\mathbf{u}^{1}=0,\quad \phi ^{1}=\phi _0^{1}(\mathbf{x}) \quad
 \text{on }\partial _{+}\Omega _1, \\
\mathbf{u}^p=0,\quad \phi ^p=0 \quad \text{on }\partial _{-}\Omega _{p}.
\end{gathered}
\end{equation}

Note that in \cite{Arjona2008} the sign of the term in $\nabla
\phi $ was the opposite. Nevertheless the techniques used there and the
present papers can be easily adapted to this stationary formulation for
proving the following result.

\begin{theorem} \label{thm6.1}
Under the same spatial regularity assumption on
the data  in Theorem 1 there exists a unique weak solution of problem
\eqref{Pinf}, \eqref{BCinf}.
\end{theorem}

\begin{proof} The existence part follows the same strategy than the proof
of part i) of Theorem \ref{thm1} and, more specifically, the proof of
\cite[Theorem 1]{Arjona2008}.
Obviously, the uncoupled problem for $\phi $
is exactly the same and the uncoupled problem for $\mathbf{u}$ is treated in
same manner (since, at this moment,  $\nabla \phi $ is assumed to be
prescribed). We recall that the coerciveness of the bilinear form $
a_{u}:V_{u}\times V_{u}\to \mathbb{R}$ given by \eqref{a_u}
requires either some assumption on the coefficients (see (69) [also denoted
as assumption $H(\rho ,\mu ,\nu )$] of \cite{Arjona2008}) or the
application of the ``dilatant argument'':  we make the changes of spatial
variables $\mathbf{y=}\lambda \mathbf{x}$, we remark that the constant in
the Poincar\'{e}'s inequality can be assumed to depend linearly on $\lambda $
(since such a constant only depends on the diameter of $\Omega $) and then,
by introducing the change of unknown
\begin{equation*}
\mathbf{v(y)}=\mathbf{u}( \mathbf{x}) =\mathbf{u}( \frac{
\mathbf{y}}{\lambda })
\end{equation*}
we see that the new condition  in terms of the new coefficients always
holds if we take  $\lambda $ large enough.

Moreover, the a priori estimates used to justify the passing to
the limit process, $n\rightarrow +\infty $, remains valid word by word since
the estimate of the right hand side of the equation for $\mathbf{u}^{n}$
only uses the norm of the vector $\| \nabla \phi ^{n-1}\| $
and thus its sign information is not relevant to this purpose. The
application of \cite[Lemma 1]{Arjona2008} (note the
resemblances with the argument on $\delta <1$ used in the proof of part (i)
of Theorem \ref{thm1} of the present paper) ends the proof of the existence of
solutions.

Concerning the proof of the uniqueness of solutions of \eqref{Pinf}, \eqref{BCinf}
we use what we call as ``cancellation method'' but now
after multiplying the equation of $\mathbf{u}^i$  by
$4\pi \rho ^i \mathbf{u}^i$ and the one of  $\phi ^i$ by
$(\rho ^i/\mu ^i)\phi ^i$ (as done in \cite{Arjona2008}).
Then adding the resultant equations, after applying Green's formula, we obtain
\begin{equation*}
\sum_{i=1}^p4\pi G\rho ^ia(\mathbf{u}^i,\mathbf{u}^i)-\sum_{i=1}^p
\frac{\rho ^i}{\mu ^i}\int_{\Omega _i}| \nabla \phi
^i( \mathbf{x}) | ^2\,d\mathbf{x}=0
\end{equation*}
(compare this with \cite[equation (52)]{Arjona2008} in which
the sign of the second term is the opposite one). Nevertheless,  by using
the standard estimate for elliptic equations (since
$\operatorname{div}\mathbf{u}^i\in H^{-1}(\Omega _i))$, and
Poincar\'{e}'s inequality we arrive to the
uniqueness conclusion if $C(\Omega _i)$ is small enough.
More exactly, the conclusion holds if $\Omega $ is such that
\begin{equation}
4\pi G\sum_{i=1}^p\rho ^i>\sum_{i=1}^p\frac{\rho ^iC(\Omega _i)K}{
\mu ^i}.  \label{critical omega}
\end{equation}

Once again, we can apply the ``dilatation argument'' in the
sense that if condition \eqref{critical omega} does not hold then we make
the dilatation $\mathbf{y=}\lambda \mathbf{x}$, the change of unknown
$\mathbf{v(y)}=\mathbf{u}( \mathbf{x}) =\mathbf{u}( \frac{
\mathbf{y}}{\lambda }) $ and we see that the last term of the right
hand side remains constant in $\lambda $ but the third term of the left hand
side appears multiplied by $\lambda $. Thus, the corresponding condition
(similar to \eqref{critical omega}) holds by taking $\lambda $ large enough
\begin{equation}
\lambda >\frac{\sum_{i=1}^p\frac{\rho ^iC(\Omega _i)K}{\mu ^i}}{
4\pi G\sum_{i=1}^p\rho ^i},  \label{dilatation}
\end{equation}
and thus the uniqueness in now proven without condition
\eqref{critical omega}, which completes the proof.
\end{proof}

 \section*{Discussion and conclusion}

A rigorous well-posedness proof of the viscoelastic-gravitational model has
been presented in this paper. The existence and uniqueness of solutions
representing a layered Earth have been carried out. For that, some techniques
of the weak solutions of partial differential equations theory have been applied.
Moreover, we have given a constructive proof of the existence of the weak solutions.
There is a clear geophysical need of this kind of models for interpretation of
observed displacement and gravity changes at volcanic areas (see introduction and,
e.g., \cite{Bonafede1998,Fernandez2001a,Folch2000}).
In previous results obtained by other authors \cite{Fernandez2001a},
introducing viscoelastic properties in all or part of the medium can extend
the effects (displacements, gravity changes, etc.) considerably and therefore
lower (and more realistic) pressure increases are required to model given
observed effects. The viscoelastic effects seem to depend mainly on the
rheological properties of the layer (zone) where the intrusion is located,
rather than on the rheology of the whole medium. Those results should be confirmed
with the described model.
Additionally, normally purely elastic half-space models are used to interpret
displacements and gravity data in active volcanic areas. Elastic-gravitational
models allow the computation of gravity, deformation, and gravitational potential
changes due to pressurized magma cavities and intruding masses together
\cite{Charco2006}, taking into account the mass interaction with the
self-gravitation of the Earth through coupling between model equations.
In \cite{Charco2006} a dimensional analysis of the elastic-gravitational model
estimating the magnitude of intrusion mass and coupling effects at the space
scale associated with volcano monitoring is performed. They show that the
intrusion mass cannot be neglected in the interpretation of gravity changes
while displacements are primarily caused by pressurization.
Therefore the intrusion of mass, together with the associated pressurization
of the magma chamber, produces distinctive changes in gravity that could be
used to interpret gravity changes without ground deformation and vice versa,
depending on what type of source plays the main role in the intrusion process.
Their theoretical experiments indicate that mass and self-gravitation could
produce changes in the magnitude and pattern of predicted gravity that may be
above microgravity accuracy and the elastic-gravitational model is a refinement
of purely elastic models which can better interpret gravity and deformation
changes in active volcanic zones. Similar studies should be done for the
viscoelastic-gravitational case described here checking the existence or
not of effects as observed in the viscoelastic-gravitational problem for
faulting (e.g., the introduction of a long-wavelength component into the time
deformation and the need of a proper consideration of gravity for near-field
computations and long time periods \cite{Fernandez1996,Fernandez2004}).
The iterative scheme presented in this work can be useful to construct a
numerical method to compute the coupled effects of gravity and viscoelastic
deformations produced by possible sources embedded in an inhomogeneous Earth.

\subsection*{Acknowledgments}

The authors thank to Professors J. Fern\'andez and J. B. Rundle for different
conversations on the modeling of the problem. The research of A. A. was partially
supported by the National Research Fund of
Luxembourg (AFR Grant 4832278).
The research of J. I. D. was partially supported
by the project ref. MTM2011- 26119 of the DGISPI (Spain) and the Research
Group MOMAT (Ref. 910480) supported by UCM. He has received also support from
the ITN FIRST of the Seventh Framework Program of the European Community
(grant agreement number 238702).


\begin{thebibliography}{00}

\bibitem{Aki-Richards} K. Aki, P.~G. Richards;
\newblock {\em Quantitative Seismology, University Science Books}.
\newblock Sausalito, California, 2002.

\bibitem{Arjona2007} A. Arjona, J.~I. D\'iaz;
\newblock Stabilization of a Hyperbolic/Elliptic system modelling the viscoelastic-gravitational deformation in a multilayered earth.
\newblock To appear in {\em Proceedings of The 10th AIMS Conference on Dynamical Systems,
Differential Equations and Applications}, Madrid, Spain, 2014.

\bibitem{Arjona2008} A. Arjona, J.~I. D\'iaz, J. Fern\'andez, J.~B. Rundle;
\newblock On the Mathematical Analysis of an
Elastic-gravitational Layered Earth Model for Magmatic Intrusion: The
Stationary Case.
\newblock {\em Pure Appl. Geophys.}, 165, 1465-1490, 2008.

\bibitem{Arjona2011} A. Arjona, T. Chac\'on, M. G\'omez;
\newblock Numerical modeling of volcanic ground deformation.
\newblock {\em(Actas de XXII
CEDYA (XI Congreso de Matem\'{a}tica Aplicada), Palma de Mallorca}, 2011.

\bibitem{Battaglia2004} M. Battaglia, P. Segall;
\newblock  The interpretation of gravity changes and crustal deformation
in active volcanic areas.
\newblock {\em Pure Appl. Geophys.}, 161, 1453-1467.

\bibitem{Battaglia2008} M. Battaglia, J. Gottsmann, D. Carbone,
J. Fern\'andez;
\newblock 4D volcano gravimetry.
\newblock {\em Geophysics}, 73, 6, WA3-WA18, 2008.

\bibitem{Battaglia2009} M. Battaglia, D. P. Hill;
\newblock Analytical modeling of gravity changes and crustal deformation at volcanoes:
The Long Valley caldera, California, case study.
\newblock {\em Tectonophysics}, 471. 45, 57. 2009.

\bibitem{Bonafede1998} M. Bonafede, M. Mazzanti;
\newblock  Modelling
gravity variations consistent with ground deformation in the Campi Flegrei
caldera (Italy).
\newblock {\em J. Volcanol. Geotherm. Res.}, 81, 137-157, 1998.

\bibitem{Brezis1984} H. Br\'ezis;
\newblock {\em Functional Analysis, Sobolev Spaces and Partial Differential Equations}.
\newblock Springer, New York,  2010.

\bibitem{Camacho2011} A. G. Camacho, P. J. Gonz\'alez, J. Fern\'andez, G. Berrino;
\newblock Simultaneous inversion of surface
deformation and gravity changes by means of extended bodies with a free
geometry: Application to deforming calderas.
\newblock {\em J. Geophys. Res.}, 116, B10401, doi: 10.1029/2010JB008165, 2011.

\bibitem{Charco2006}M. Charco, J. Fern\'andez, F. Luz\'on, J. B. Rundle;
\newblock On the relative importance of self-gravitation and elasticity
in modeling volcanic ground deformation and gravity changes.
\newblock {\em Journal of Geophysical Research},
111, B03404, doi: 10.1029/2005JB003754, 2006.

\bibitem{Charco2014} M. Charco, P. Gal\'{a}n del Sastre;
\newblock Efficient inversion of 3D Finite Element models of volcano
deformation.
\newblock {\em Geophys. J. Int.}, 196, 1441-1454, doi:10.1093/gji/ggt490, 2014.

\bibitem{Diaz2002} J.~I. D\'{\i}az, R. Quintanilla;
\newblock  Spatial and continuous dependence estimates in linear viscoelasticity.
\newblock {\em J. Math. Anal. And Applic.}, 273, 1-16, 2002.

\bibitem{Duvaut1972} G. Duvaut, J. L. Lions;
\newblock {\em Les in\'equations en m\'ecanique et en physique}.
\newblock Dunod, Paris, 1972.

\bibitem{Fernandez1994a} J. Fern\'andez, J. B. Rundle;
\newblock Gravity changes and deformation due to a magmatic intrusion in a two-layered
crustal model.
\newblock {\em J.Geophys.Res.}, 99, 2737-2746, 1994a.

\bibitem{Fernandez1994b}  J. Fern\'andez, J.~B. Rundle;
\newblock FORTRAN program to compute displacement, potential and gravity changes due
to a magma intrusion in a multilayered Earth models.
\newblock {\em Comput.Geosci.} , 23, 231-249, 1994b.

\bibitem{Fernandez1996} J. Fern\'andez,  Ting-To Yu, J.~B. Rundle;
\newblock Horizontal viscoelastic-gravitational displacement due to a
rectangular dipping thrust fault in layered Earth model.
\newblock {\em J. Geophys.Res.},
101, 6, 13.581-13.594. (correction: Journal of Geophysical Research, 103,
B12, 30283-30286, 1996.

\bibitem{Fernandez1999} J. Fern\'andez, J. M. Carrasco, J. B. Rundle, V. Ara\~na;
\newblock Geodetic methods for detecting volcanic unrest:
a theoretical approach.
\newblock {\em Bull. Volcanol.} , 60, 534-544, 1999.

\bibitem{Fernandez2001a} J. Fern\'andez, K. F. Tiampo, J. B. Rundle;
\newblock Viscoelastic displacement and gravity changes due to point
magmatic intrusions in a gravitational layered solid earth.
\newblock  {\em Geophys. J.
Int.}, 146, 155-170, 2001a.

\bibitem{Fernandez2001b} J. Fern\'andez, M. Charco, K.~F. Tiampo, G. Jentzsch,
 J.~B. Rundle;
\newblock Joint interpretation of
displacement and gravity data in volcanic areas. A test example: Long Valley
Caldera, California.
\newblock  {\em J. Volcanol.Geotherm. Res.}, 28, 1063-1066, 2001b.

\bibitem{Fernandez2001c} J. Fern\'andez, K.~F. Tiampo, G. Jentzsch, M. Charco,
 J. B. Rundle;
\newblock Inflaction or deflation?. New
results for Mayon volcano applying elastic-gravitational modeling.
\newblock  {\em Geophys.Res.Lett., 28, 12, 2349-2352, 2001c}.

\bibitem{Fernandez2004} J. Fern\'andez, J. B. Rundle;
\newblock Postseismic visoelastic-gravitational half space computations:
Problems and solutions.
\newblock {\em Geophys. Res. Lett.}, 31, 2004.

\bibitem{Folch2000} A. Folch, J. Fern\'andez, J.~B. Rundle, J. Mart\'{\i};
\newblock Ground deformation in a viscoelastic medium composed
of a layer overlying a half space: A comparison between point and extended
sources.
\newblock {\em Geophys.J.Int.}, 140, 37-50, 2000.

\bibitem{Fung1965} Y.~C. Fung;
\newblock {\em Foundations of Solid Mechanics.}
 \newblock Prentice-Hall, Englewood Cliffs, N.J., 1965.


\bibitem{Gilbarg1977} D. Gilbarg, N.~S. Trudinger;
\newblock {\em Elliptic Partial Differential Equations of Second Order.}
\newblock Springer-Verlag, Berlin, 1977.

\bibitem{McTigue1987} D.~F. McTigue;
\newblock Elastic stress and
deformation near a finite spherical magma body: Resolution of the point
source paradox.
\newblock {\em J. Geophys.Res.}, 92, 12931-12940, 1987.

\bibitem{Mogi1958} K. Mogi;
\newblock Relations between the
eruptions of various volcanoes and the deformation of the ground surfaces
around them.
\newblock {\em Bull.Earth.Res.}, 92, 12 931-12 940, 1958.

\bibitem{Okada1958} Y. Okada;
\newblock Surface deformation due
to shear and tensile faults in a half-space.
\newblock {\em Bull. Seism. Soc. Am.}, 75, 1135-1154, 1958.

\bibitem{Okubo1992} S. Okubo;
\newblock Gravity and potential
changes due to shear and tensile faults in a half-space.
\newblock {\em J. Geophys.Res.}, 97, 7117-7144, 1992.

\bibitem{Pollitz1997} F.~F. Pollitz;
\newblock Gravitational
viscoelastic postseismic relaxation in a layered spherical earth.
\newblock {\em J.
Geophys.Res.}, 102, 17921-17941, 1997.

\bibitem{Rundle1980} J. B. Rundle;
\newblock Static
elastic-gravitational deformation of a layared half space by point couple
sources. \newblock {\em J. Geophys.Res.}, 85, 5355-5363, 1980.

\bibitem{Rundle1981a} J. B. Rundle;
\newblock  Numerical
Evaluation of static elastic-gravitational deformation of a layared half
space by point couple sources. \newblock {\em Rep., Sand}, 80-2048, 1981a.

\bibitem{Rundle1981b} J. B. Rundle;
\newblock   Deformation gravity and potential changes due to volcanic loading of the crust.
\newblock {\em J.
Geophys.Res.}, 87, 10, 729-10, 744, 1981b.

\bibitem{Rundle1982b} J. B. Rundle;
\newblock  Viscoeslastic-Gravitational Deformation by a Rectangular Thrust Fault in a
Layered Earth.
 \newblock {\em J. Geophys.Res.}, 87, 9, 7787-7796, 1981b.

\bibitem{Rundle1983} J. B. Rundle;
\newblock  Correction to
"Deformation, gravity and potential changes due to volcanic loading of the
crust".
 \newblock {\em J. Geophys.Res.}, 88, 10, 647-10, 653, 1983.

\bibitem{Yang1998} X. Yang, P. M. Davis, J. H. Dieterich;
\newblock  Deformation from inflation of a dipping finite prolate spheroid in
an elastic half-space as a model for volcanic stressing.
 \newblock {\em J.Geophys.Res.},
93, 4249-4257, 1998.

\bibitem{Yu1996} T. T. Yu, J. B. Rundle, J. Fern\'andez;
\newblock  Surface deformation due to strike-slip fault in an elastic
gravitational layer overlying a viscoelastic gravitational half-space.
\newblock {\em J. Geophys. Res.}, 101, 3199-3214.
(Correction: Journal of Geophysical Research, 104, 15313-15316.), 1996.


\end{thebibliography}


\end{document}