\documentclass[a4paper,11pt]{article}
%\documentclass[dvips]{article}
%\documentclass[12pt,letterpaper]{article}
%\usepackage{times}
\usepackage{graphicx}
\usepackage{amsmath}
\usepackage{graphicx}
\usepackage{epstopdf}
\usepackage{setspace}
\usepackage[normalsize]{caption}

\usepackage[numbers,round,colon,authoryear]{natbib}
%\usepackage{natbib}


\newcommand{\ud}{{\rm d}}

\usepackage{url}


\setlength{\textwidth}{17.5cm}
\setlength{\textheight}{26.8cm}
\setlength{\hoffset}{-2.5cm}
\setlength{\voffset}{-2.9cm}


%\usepackage{pause}
%\definecolor{webgreen}{rgb}{0, 0.15, 0} % less intense green
%\definecolor{webblue}{rgb}{0, 0, 0.15} % less intense blue
%\definecolor{webred}{rgb}{0.15, 0, 0} % less intense red
%\usepackage[pdfauthor={Boris Jeremic}, 
%            colorlinks=true,
%            linkcolor=webblue,
%            citecolor=webred,
%            urlcolor=webgreen,
%            linktocpage,
%            dvips]{hyperref}
%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ovo je predefinisanje datuma za cirilicu ( vidi Kopka93 strana 175 )
\def\today
{\number\day. \space \ifcase\month\or
January\or
February\or
March\or
April\or
May\or
June\or
July\or
August\or
September\or
October\or
November\or
December\fi,\space
\number\year}

%time of day as found in layout.tex
% ----------------------------------------------------------------------
%
% TIME OF DAY
%
\newcount\hh
\newcount\mm
\mm=\time
\hh=\time
\divide\hh by 60
\divide\mm by 60
\multiply\mm by 60
\mm=-\mm
\advance\mm by \time
\def\hhmm{\number\hh:\ifnum\mm<10{}0\fi\number\mm}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



%\def\baselinestretch{2}
\def\baselinestretch{1.1}


%Mahdi
%% The lineno packages adds line numbers. Start line numbering with
%% \begin{linenumbers}, end it with \end{linenumbers}. Or switch it on
%% for the whole article with \linenumbers after \end{frontmatter}.
%\usepackage{lineno}
%\usepackage[left,pagewise,displaymath, mathlines]{lineno}
\usepackage[left,displaymath, mathlines]{lineno}

\begin{document}



\begin{center}
%{\Large \bf Simulation of Fully Coupled Behavior of Porous Media: Verification}
{\Large \bf
Solution Verification Procedures for 
Modeling and Simulation of Fully Coupled Porous Media: Static and Dynamic  Behavior}
\end{center}
\vspace*{1cm}

\begin{center}
{\large \bf Panagiota Tasiopoulou} \\
{PhD Student, 
National Technical University of Athens, 
Athens, Greece}\\ ~ \\
%Mahdi Taiebat \scalebox{1}{\includegraphics[viewport=0 4.5 0 0]{mahdi.eps}}
{\large \bf Mahdi Taiebat} \\
Associate Professor, 
Department of Civil Engineering, The University of British Columbia, 
Vancouver, BC, Canada\\       ~ \\
%%%%%%%%
{\large \bf Nima Tafazzoli} \\
Geotechnical Engineer, Tetra Tech EBA, Vancouver, BC, Canada \\ ~ \\ 
{\large \bf Boris Jeremi{\'c}}\footnote{Corresponding Author: Boris Jeremi{\'c},  \\ 
Department of Civil
and    Environmental    Engineering   University  of  California,  Davis,  CA,
U.S.A. \texttt{jeremic@ucdavis.edu}, and \\
Earth   Science  Division,  Lawrence Berkeley National Laboratory, Berkeley, CA,
U.S.A. \texttt{bjeremic@lbl.gov}}
\\
Professor, 
Department  of  Civil  and  Environmental  Engineering,  
University  of  California, Davis, CA, U.S.A. \\ and \\
Faculty Scientist, 
Earth   Science  Division,  
Lawrence Berkeley National Laboratory, Berkeley, CA, U.S.A.
\end{center}
%%%%%%%%
%and
%%%%%%%%%
%Yiorgos Perikleous\affilnum{4}

%\vspace*{2cm}
%\begin{center}
%{\sc Version: \today{} , \hhmm{} }
%\end{center}
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\newpage
%\tableofcontents
%\newpage


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%\end{abstract}

%\linenumbers

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{Abstract}
%{\sc \Large I have updated it this}

Numerical prediction of dynamic behavior of fully coupled saturated porous media
is  of  great  importance in many engineering problems. Specifically, static and
dynamic  response of soils -- porous media with pores filled with fluid, such as
air, water, etc. -- can only be modeled properly using fully coupled approaches.
Modeling  and  simulation  of  static  and  dynamic  behavior  of  soils require
significant  Verification  and  Validation  (V\&V)  procedures in order to build
credibility  and  increase  confidence  in  numerical  results.  By  definition,
Verification  is  essentially  a mathematics issue and it provides evidence that
the model is solved correctly, while Validation, being a physics issue, provides
evidence that the right model is solved.
%
This  paper  focuses  on  Verification  procedure for fully coupled modeling and
simulation  of  porous  media. Therefore, a complete Solution Verification suite
has  been  developed  consisting  of  analytical  solutions  for both static and
dynamic problems of porous media, in time domain. Verification for fully coupled
modeling and simulation of porous media has been performed through comparison of
the  numerical  solutions  with  the analytical ones. Modeling and simulation is
based on the so called, $u$--$p$--$U$ formulation.
%\citep{Zienkiewicz84, Jeremic2004d,
%Jeremic2007e}. 
Of    particular interest are numerical  dispersion  effects which determine the
level of numerical accuracy.
%
These effects   are 
investigated in detail, in an effort to  suggest a compromise between 
numerical error and computational cost.
%
%In    this    paper,    a  solution  verification   suite (a subset of the whole
%verification  suite,  divided  into  code  and solution verification processes),
%consisting  of  both static and dynamic problems, is developed for fully coupled
%modeling and simulation of porous media. Modeling and simulation is based on the
%so  called,  $u$--$p$--$U$  formulation  \citep{Zienkiewicz84, Jeremic2004d,
%Jeremic2007e}.
%
%Static  and  dynamic  examples  are  used    to  verify  the numerical solutions
%through  comparison  with analytical solutions. 
%
%Of  particular interest are numerical  dispersion  effects. 
%
%These effects   are 
%investigated in detail in an effort to  suggest a compromise between 
%numerical error and the computational cost.



\vspace*{0.3cm}

\noindent{\bf    keywords:}    verification,  finite    elements,  fully coupled
analysis, porous media, numerical dispersion
%elastic-plastic
%Note  (mahdi):  we  are not dealing with elastic plastic materials here in this
%paper -- our focus is rather only on elastic material.
\newpage








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



One    of   the  main  questions  raised  when  results  of a numerical analysis
are
evaluated      is:      ``How    much    can    (should)    we   trust    model 
implementations
and how much can (should) we trust numerical simulations?"
Trust  in  numerical  implementation is  gained through a verification process,
while,  trust in results obtained from numerical simulation, is gained through a
validation process \citep{Oberkampf2002, Roy2011}.
%


\begin{itemize}

\item[]  {\bf  Verification}  is    a  process  of determining  that    a  model
implementation
accurately  represents the developer's conceptual description and specification.
It is essentially a mathematics issue and it provides evidence that the model is
solved correctly.
%

\item[]  {\bf  Validation},  on  the other hand, is a process of determining the
degree to which
a model is accurate representation of the real world from the perspective of the
intended uses of the model. It is a physics issue, and it provides evidence that
the correct model is solved.
%

\end{itemize}

\noindent
Verification and Validation (V\&V) are  the  primary means of assessing accuracy
in  modeling and computational simulations in order to build confidence and
credibility  in  numerical  predictions.
\citet{Oberkampf2002}
mainly  focus on trying to accurately model the ``real world", aiming to provide
verified and validated
computer simulation of the reality.
%
Slightly    different  approach  is taken by \citet{Oden2010a, Oden2010b}, where
the
actual  purpose    of   numerical  simulation  is  to facilitate the engineering
decision
process,    rather  than  to explicitly represent the ``real world". Despite the
alternative  approaches on the objectives of numerical modeling, it is important
to note that both schools of thought
place significant emphasis on the need
for
validation of numerical predictions under uncertainty. In an attempt to describe
how  V\&V  provide  the  connection  between  reality  and  numerical simulation
(computer  implementation) in the process of reaching an engineering decision, a
flow chart has been constructed and shown in Figure \ref{RoleVV}. 
%
\begin{figure}[!htb]
\begin{center}
{\includegraphics[width=7.0cm]{FIGS/RoleVV_NEW_knowledge.pdf}}
%\vspace*{-0.5cm}
\caption{\label{RoleVV} Schematic representation of the role of Verification and
Validation  on  numerical  modeling  (inspired     from     \protect{\citet{Oberkampf2002,Roy2011}}     and
\protect{\citet{Oden2010a}}).}
\end{center}
\end{figure}
%




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%





%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%
%\paragraph{Verification Procedures.}
%\label{VandV_P}


It  is  evident  that  Verification  and Validation procedures have separate but
complementary  objectives  by  definition:  i)  the  model  is  solved correctly
(Verification)  and ii) the correct model is solved (Validation). On this basis,
they  should  be  dealt as two distinct stages requiring different activities in
order   to   achieve  them.  In  particular,  two  main  activities  related  to
verification procedure can be described as \citep{Oberkampf2002}:

\begin{itemize}


         \item  {\bf Code  Verification} which identifies  and  removes errors in computer
         coding.  It  can  further  be  separated into:
%
         (a)  {\bf  numerical  algorithm  verification},  used  to  ensure  that all
         the numerical
         algorithms are correctly implemented in the code, and
%
         (b)  {\bf software  quality  assurance},  used to ensure that the program
         system is implemented reliably so that it produces repeatable
         (stable)  results  with  different  compilers  and  on  different  computer
         architectures.



   \item  {\bf  Solution Verification} which quantifies  numerical errors in
   computed  solutions.  Numerical accuracy can only be estimated by direct comparison
   with analytical solutions for a given application. Thus, evidence that the model can be solved correctly is provided.



     \end{itemize}
  
On  the  other  hand,  Validation  procedure  involves  numerical  simulation of
realistic  conditions  and  thus, direct comparison with experiments and/or well
documented  case  histories  is  needed.  It  is  common  in  practice that only
validation  procedures  are  performed,  considered  misleadingly to show enough
evidence  that  the  numerical  code is appropriate to solve correctly realistic
numerical  models, while the verification stage is omitted. Therefore, it should
be  emphasized  that  verification  is the initial necessary step to ensure that
numerical  solution of the model is accurate while validation is the final stage
to  ascertain  the  numerical model can represent realistic conditions and thus,
numerical predictions can be trusted.


The  current  study  focuses  on the Verification stage and particularly, on the
solution  verification  procedure,  assuming  that code verification, which is a
strictly programming issue, has been successfully performed. The twofold goal of
this  paper  is  to: a) provide a complete solution verification suite for fully
coupled,  static  and dynamic behavior of porous solid-fluid systems, consisting
of   analytical  solutions  for  porous  media  in  time  domain,  available  in
literature,  which  can  be  used to verify any finite element/finite difference
code  and  b)  to  use  this  particular solution verification suite in order to
verify   the   UCD  Computational  Geomechanics  Group's  implementation  of  3D
$u$--$p$--$U$   finite  elements  \citep{Jeremic2007e},  which  is
available through an open source license. The available analytical solutions are
provided  for  idealized  problems,  based on assumptions such as linear elastic
behavior  and  simplified geometries (mostly 1D). Validation of coupled behavior
using  elasto-plastic  material  models,  is  the  subject  of another paper
\citet{Tasiopoulou2014b}.



   It  is  important  to  note  that  this paper provides, to our knowledge, the
   first verification of fully coupled (porous solid -- pore fluid) modeling and
   simulation.  In  the  past, there were few published papers that used some of
   the  analytical  and numerical solutions, used here, for comparison purposes.
   For  example  a  paper by \citet{Gajo1994}, published a while ago, focused on
   comparing  different  methods,  namely u-p-U, u-U and u-w and used two closed
   form  solutions  for  that  purpose. However, they did not provide details of
   numerical  algorithms used, since they used the same algorithm and parameters
   (whatever they were) comparing different methods (as noted above, u-p-U, u-U
   and  u-w).  Focus  of this paper is entirely on verification of u-p-U method.
   For  this  purpose all the available closed form solutions (to our knowledge)
   are    used  and  all  the   algorithms and parameters that we used are fully
   documented. It can
   be  claimed  that  the  u-p-U formulation developed by \citet{Zienkiewicz84},
   together with our implementation (presented verification, and validation provided by
   \citep{Tasiopoulou2014b}), is now complete.





%It      should    be  mentioned  that the current verification  process has been
%performed using the deterministic approach. While this might be appropriate
%for  solution  verification  and  code verification, validation must address the
%issue   of  uncertainty,  as  noted  by  \citet{Oden2010a, Oden2010b, Oden2011}.
%Such  undertaking  is  planned  for  near future and will rely on Probabilistic
%Elastic-Plasticity and Stochastic Elastic-Plastic Finite Element Method
%\citep{Jeremic2005a, Sett2005a, Sett2009b, Sett2010a}.





%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%
% \clearpage






%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%

%\section{Theoretical-Numerical  Formulation  
%%of  Porous Media 
%and Implementation
%Background}
%\label{Formulation}


%{\sc \Large THIS NEEDS TO BE UPDATED, it is from our old paper... The following
%subsection  has  been  changed  and  the  constitutive integration part has been
%omitted..}


%\subsection{Governing Equations of Porous Media}
\section{Governing Equations of Porous Media}
\label{Governing Equations of Porous Media}


The response of fully saturated soils under transient loads plays  an  important
role  in  geotechnical  engineering.  A  literature  review  on  the  subject of
two-phase  materials  identifies  two  major  theories which have been developed
during  the last century and are currently used, namely the Biot theory  and
the      theory      of      Porous      Media .      Biot's   theory 
%\citep{Biot1941,Biot1956a,Biot1956b,Biot1956c} 
is initially based on
the    work    of    \citet{Terzaghi23},    who    proposed    a    model   for 
one-dimensional
consolidation.  Theory      of      Porous      Media  is  based on the axioms of continuum theory of
mixtures      extended      by      the      concept      of    fractions    by 
\citet{Bowen1980,Bowen1982,Ehlers1993}.

\citet{Biot1941}  presented  a theoretical description of linear elastic porous
materials  fully  saturated  by  a  viscous fluid under quasi-static conditions.
Extensions  of  his  theory  to  anisotropic  cases  \citep{Biot1955}  and  to
poroviscoelasticity  \citep{Biot1956a,Biot1956b,Biot1956c}  followed. In Biot's
theory  a  fully  saturated  material  is  considered,  which  consists  of  two
compressible  parts:  the  solid  skeleton  and  the pore fluid. Thus, the solid
displacements,  $u_i$  and  the  fluid displacements, $U_i$, are introduced. The
equations  of  motion are formed by stating the balance of momentum, or else the
dynamic  equilibrium, for the solid skeleton and the pore fluid separately.
\citet{Zienkiewicz84}  presented  a  slightly  different formulation,
which  is  based on the concept of Biot's theory; however, it is more convenient
for  numerical  solution.  Alternatively  to the fluid displacements, a relative
displacement  of  fluid to the solid, $w_i$ (the so-called seepage displacement)
is  used. Moreover, no additional mass density, $\rho_\alpha$, is introduced for
the  dynamic  interaction  between  the  fluid  and the skeleton, in contrast to
Biot's  theory. 

Formulation   presented by  \citet{Zienkiewicz84} is the most general one and is
used in our developments.
%
They did not, however, provide verification of
developed  formulation    that  will  prove  that the formulation is implemented
accurately  and  that  the examples are solved correctly. 
%
Detailed  solution  verification of $u-p-U$ formulation is precisely the goal of
this paper.



Fully coupled ($u-p-U$) formulation is given below in
some  detail,  and is needed as reference for some of the verification solutions
given later in the paper.
%
We  start from the  relationship  between  effective  stress,
total  stress and pore pressure, written as $\sigma^{'}_{ij}= \sigma_{ij}+\alpha
\delta_{ij}      p$,      where   $\sigma^{'}_{ij}$  is  the  effective  stress 
tensor\footnote{stresses are positive in tension}, $\sigma_{ij}$ is total stress
tensor,  $\delta_{ij}$  is  Kronecker  delta, $p$ is the pore fluid pressure\footnote{pore pressure is positive in compression} and
$\alpha$  is  the  Biot constant that depends on the geometry of material voids.
For  the  most  part,  in  soil  mechanics  problems,  $\alpha \approx 1$ can be
assumed,  so  that  the  relationship between total and effective stress becomes
$\sigma^{'}_{ij}   =   \sigma_{ij}+\delta_{ij}  p$,  which  corresponds  to  the
classical effective stress definition by \citet{Terzaghi43}.
                                                                                                                                                                                                                                 
The    material  behavior of soil  skeleton (porous solid) is  fully  dependent 
on the pore fluid pressures. The behavior of pore fluid is assumed to be elastic
and thus all the material nonlinearity is concentrated on the soil skeleton. The
soil  behavior ('mixture' of soil skeleton and pore fluid) can thus be described
using  single--phase  constitutive  analysis approach for skeleton combined with
the  full  coupling  with pore fluid. The overall equilibrium (momentum balance)
equation for the soil-fluid ``mixture" can be written as
%%%
\begin{equation}
   \sigma_{ij,j} - \rho \ddot{u}_i - \rho_f \ddot{w}_i + \rho b_i = 0
\label{upU_GEequilibrium}
\end{equation}
%%%
where  $\ddot{u}_i$  is  the  acceleration  of the solid part, $b_i$ is the body
force  per  unit  mass,  $\ddot{w}_i$  is the fluid acceleration relative to the
porous solid.
For fully saturated porous media (no air trapped inside),
density is equal to $\rho = n \rho_f + (1 - n) \rho_s$,
where $n$ is the porosity,
$\rho_s$ and $\rho_f$ are the soil particle and fluid density, respectively.
For the pore fluid, the equation of momentum balance can be written as
%%%
\begin{eqnarray}
 -p_{,i} - {R_i} - \rho_f \ddot{u}_i - \frac{\rho_f \ddot{w}_i}{n} +\rho_f b_i =
 0
\label{GEequilibriumFluid}
\end{eqnarray}
%%%
where $R$ is the viscous drag force. According to the Darcy's seepage law,
the viscous drag forces $R$ between soil matrix and pore fluid can be written as
${R_i} = k_{ij}^{-1}  \dot{w}_j$,
where
$k_{ij}$  is  the   tensor of anisotropic intrinsic  permeability  coefficients.
For
simple case of isotropic permeability, scalar value of permeability $k$ is used,
so as ${R_i} = k^{-1} \dot{w}_j$.
The        permeability   ${k}$ used  here with dimension of [$L^{3}TM^{-1}$] is
different from Darcy permeability (hydraulic conductivity), ($k_{D}$), which has
the dimension of velocity, i.e. [$LT^{-1}$].
Their values are related by $k = k_{D} / g \rho_f$,
where  $g$  is  the  gravitational  acceleration  and $\rho_f$ is the density of
the pore fluid.



The final equation is the mass conservation of the fluid flow expressed by
%%
\begin{eqnarray}
 \dot{w}_{i,i} + \alpha \dot{\varepsilon}_{ii} + \frac{\dot{p}}{Q} = 0
\label{GEconservationFlow}
\end{eqnarray}
%%
where bulk stiffness  of the mixture $Q$ is expressed as
${1}/{Q} = {n}/{K_f} + {(\alpha-n)}/{K_s}$,
and $K_s$ and $K_f$ are the bulk moduli of the solid (particles, not the skeleton) and fluid phases, respectively.

In    the    above    governing  equations,  convective terms of lower order are
omitted. A change of variable is performed by introducing an
alternative variable $U_i$, defined as $ U_i = u_i + {w_i}/{n}$, that represents
absolute  displacement  of  the  pore  fluid.  Thus,  the  basic set of unknowns
comprises  the  soil  skeleton displacements $u_i$, the fluid pore pressure $p$,
and  the  fluid displacements $U_i$. The set of the modified governing equations
is summarized below:
%
\begin{eqnarray}
\sigma_{ij,j}^{'}-(\alpha-n) p_{,i}+(1-n) \rho_s b_i-(1-n) \rho_s \ddot{u}_i +
n R_i=0
\label{34}
\end{eqnarray}
\begin{eqnarray}
-n p_{,i}+n \rho_f b_i-n \rho_f \ddot{U}_i
- n R_i=0
\label{35}
\end{eqnarray}
\begin{eqnarray}
-n \dot{U}_{i,i}=(\alpha-n) \dot{\varepsilon}_{ii}+\frac{1}{Q} \dot{p}
\label{36}
\end{eqnarray}

From the modified equation set (\ref{34}), (\ref{35}), and (\ref{36}), we can see that only
$\ddot{u}_i$  occurs in the first equation, and only $\ddot{U}_i$ in the second,
thus leading to a convenient diagonal form in discretization.

This  theoretical  approach describes the interaction of solid skeleton and pore
fluid  with  three  major  equations,  governing  both dynamic and quasi--static
phenomena:     i)    equilibrium    of    total    solid-fluid    ``mixture"    -- 
Eqs.~(\ref{upU_GEequilibrium}) and  (\ref{34}), ii) equilibrium of pore fluid -- 
Eqs.~(\ref{GEequilibriumFluid}) and  (\ref{35}),  and  iii) mass conservation of
the  fluid  flow  -- Eqs.~(\ref{GEconservationFlow}) and (\ref{36}). The boundary
conditions imposed on these variables will complete the problem. These equations
together  with an appropriate constitutive relationship describe the behavior of
solid  skeleton  accounting  for  its full interaction with the pore fluid under
both  quasi-static and dynamic conditions. Herein, the constitutive relationship
between  increments of effective stress and strain of the solid skeleton is given in a general
incremental form:
%
\begin{equation}
  d\sigma^{'}_{ij} = D_{ijkl} d\varepsilon_{kl}
\label{constitutive}             
\end{equation}
%%%
where $\varepsilon_{ij} = (u_{i,j}+u_{j,i})/2$ is the small strain tensor of the
solid  skeleton,  and $D_{ijkl}$ is the tangent stiffness that can be elastic or
elasto-plastic.

%In order to solve the equation of motion, the degrees of freedom must  be  determined. In retrospect, nowadays three general continuum formulations \citep{Zienkiewicz84} are common for modeling of the fully coupled problem %(soil skeleton -- pore fluid) in  geomechanics, namely the:
%
%\begin{enumerate}
%\item A combination of the pore pressure $p$ and the solid displacement $u_i$ with four unknowns in 3D, called $u$--$p$ (two--field formulation).
%\item The solid displacement $u_i$ and the fluid displacement $U_i$ with six unknowns in 3D, called $u-U$ (two--field formulation).
%\item A combination of the pore pressure $p$, the solid displacement $u_i$ and the fluid displacement $U_i$ with seven unknowns in 3D, called $u$--$p$--$U$ (three--field formulation).
%\end{enumerate}
%
%
%The $u$--$p$ formulation captures the movements of the soil skeleton and the change of the pore pressure. Despite the fact that it is the most simplistic one of the three mentioned above, it can be valuable when the movements of %the water are not important. This formulation neglects the relative accelerations of the pore fluid (except for combined (same) acceleration of pore fluid and solid), and in one version neglects the compressibility of  the %fluid  (assuming complete incompressibility of the pore fluid). This formulation must rely on Rayleigh damping to model velocity proportional energy dissipation (damping). The  $u-U$  formulations  tracks the movements of both %the soil skeleton and the pore  fluid.  This  formulation  is complete in the sense of basic variable, but might still experience numerical problems (volumetric locking) if the difference in volumetric compressibility of the %pore fluid and the solid skeleton is large.

The  $u$--$p$--$U$  formulation  considers  compressible  pore fluid and solid grains.
Additionally,  due  to  computation of fluid inertia forces, this formulation is
applicable  to  any  range  of  frequency  for  which the hypothesis of constant
permeability  is  valid  \citep{Gajo1994}.  Moreover,  the $u$--$p$--$U$ formulation
resolves the issues of volumetric locking by including the displacements of both
the solid skeleton and the pore fluid, and the pore fluid pressure as well. This
formulation  uses (dependent) unknown field of pore fluid pressures to stabilize
the  solution  of  the coupled system. The pore fluid pressures are connected to
(dependent   on)   displacements  of  pore  fluid,  as,  with  known  volumetric
compressibility of the pore fluid, pressure can be calculated.

An    important    advantage    of $u$--$p$--$U$  formulation over commonly used
$u$--$p$   formulations   is  that  velocity  proportional  viscous  damping  is
introduced  directly  through the damping tensor which is a function of porosity
and  permeability  of  the  soil  skeleton.  This  viscous  damping provides for
physically  based, velocity proportional energy  dissipation,  stemming from the
interaction  of  pore fluid and the porous solid (soil) skeleton. In addition to
these   advantages,  the  inclusion  of  both  solid  skeleton  and  pore  fluid
displacements  in  the  field  of  unknowns  allows for independent treatment of
accelerations  of both constituents (skeleton and fluid) which improves accuracy
of  simulations.  
%
Some   details   of   the   finite   element   formulation   and  implementation
is  given  in  the  appendix for reference.





%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%\clearpage
\section{Solution Verification Suite}
\label{Verification_Suite}


Solution     verification    of    the  $u$--$p$--$U$    formulation    and  its
implementation  is  presented  in the following. Analytical solutions for porous
media,  found  in literature, are used to verify numerical solutions obtained by
application  of  the  developed  code.  The analytical solutions in time domain,
available in literature, have been developed for linear elastic porous media and
simplified geometries (mostly 1D) and they are divided into two main categories,
based  on:  a)  consolidation  theory,  where  no inertia effects are taken into
account,  called  quasi-static  poroelasticity,  and  b)  dynamic version of the
two-phase theory, called poroelastodynamics. Characteristics solutions from both
categories  are  chosen  in  the  following  to  verify the $u$--$p$--$U$ finite
element formulation, as developed by \cite{Jeremic2007e}.

It  should  be  mentioned  that  the  majority  of  the analytical solutions are
provided  for  idealized  examples,  such as 1D linear elastic wave propagation,
which  are  not always representative of realistic cases in terms of 3D geometry
and  non--linearity.  However,  while  they  would  seem  to  be  inadequate for
Validation  purposes,  which require realistic experiments and/or well documented
case      histories      (see \citet{Tasiopoulou2014b} for more details), they successfully
serve  the  mathematical objective of Verification, ensuring the accuracy of the
numerical solution provided by the code.
                                                       


\subsection{Quasi--static Poroelasticity}  
\label{Quasi--static Poroelasticity}

This  section  focuses on analytical solutions based on coupled theory of linear
porous   media   under   quasi--static  conditions  \citep{Coussy95,Coussy2004}.
Eqs.~(\ref{upU_GEequilibrium}) and (\ref{GEequilibriumFluid}), which are based
on  equilibrium of total ``mixture" and fluid respectively, are reformed in order
to  apply  for  quasi--static  phenomena  (inertia  terms  are  omitted).  Thus,
Eq.~(\ref{upU_GEequilibrium}) is written as:
\begin{equation}
  \sigma_{ij,j} = 0 
\label{equil_mix_quasi}             
\end{equation}
and Eq.~(\ref{GEequilibriumFluid}), assuming isotropic permeability, is written as:
\begin{eqnarray}  
 -p_{,i} - \frac{\dot{w}_i}{k} +\rho_f b_i = 0  
\label{equil_Fluid_quasi}                
\end{eqnarray}
Differentiating Eq.~(\ref{equil_Fluid_quasi}) with respect to spatial coordinates, we obtain:
\begin{eqnarray}  
 \dot{w}_{i,i} = -k p_{,ii}
\label{equil_Fluid_quasi_darcy}                
\end{eqnarray}
Combining              Eqs.~(\ref{equil_Fluid_quasi_darcy})              and
(\ref{GEconservationFlow}), the mass conservation equation becomes:
\begin{eqnarray}  
 \alpha \dot{\varepsilon}_{ii} + \frac{\dot{p}}{Q} = k p_{,ii}
\label{mass_cons_quasi}                
\end{eqnarray}
A simple linear constitutive relationship between effective stress and strain is
considered,
\begin{eqnarray}  
 \sigma_{ij} = (K - \frac{2}{3}\mu)\varepsilon_{v}\delta_{ij}+2\mu\varepsilon_{ij}+\alpha p\delta_{ij}
\label{shear_strain}                
\end{eqnarray}
where  $K$ is the bulk modulus, $\mu$ is the shear modulus, and $\varepsilon_{v}$ is
the    volumetric    strain   of   porous   solid   skeleton.   Combination   of
Eqs.~(\ref{shear_strain})   and   (\ref{equil_mix_quasi}),  and   subsequent
differentiation with respect to spatial coordinates, leads to $(K + (4/3)\mu)\varepsilon_{v,ii} -
\alpha p_{,ii} = 0$, which can be solved for $\varepsilon_{v}$ as:
\begin{eqnarray}  
 \varepsilon_{v} = \frac{\alpha p}{K + (4/3)\mu}
\label{strain_p}                
\end{eqnarray}

Lastly, differentiating $\varepsilon_{v}$ from Eq.~(\ref{strain_p}) with respect to time and substituting it in Eq.~(\ref{mass_cons_quasi}) in combination with the relationship $K_{u}=K +\alpha^{2}Q$ (where $K_u$ is the undrained bulk modulus of soil), the diffusion equation is derived as:
\begin{eqnarray}  
 \dot{p} = c_{f} p_{,ii}
\label{diffusion}                
\end{eqnarray}
where $c_f$ is the fluid diffusivity coefficient, defined as:
\begin{eqnarray}  
 c_{f} = kQ\frac{K+(4/3)\mu}{K_{u}+(4/3)\mu}
\label{c_f}                
\end{eqnarray}

\citet{Coussy95,Coussy2004} has presented analytical solutions for the diffusion
Eq.~(\ref{diffusion})  for different boundary condition problems which have been
used as verification examples and are studied in the following.



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\subsubsection{Vertical Consolidation of a Soil Layer}
\label{Vertical_Consolidation_of_a_Soil_Layer}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
{\bf Analytical Solution by \citet{Coussy2004}} -- 
%
%
%he    consolidation   process  can  be  defined  as  follows:  When  an elastic
%soil layer is subjected to an external change in mean normal stress, immediately
%the  water will alone sustain this increment of mean normal stress and cause the
%build-up  the  excessive pore water pressure. In the progress of the flow of the
%water to the surface, the load is gradually transferred to the soil skeleton and
%the  excessive  pore  water  pressure  will  dissipate.  At  the  same time, the
%settlement of the soil layer occurs. As settlement is usually a major concern in
%geotechnical engineering, this is a key problem in soil mechanics.
%
A   soil  layer of thickness  $h$ in the $z$ direction, composed  of  isotropic,
homogeneous   and  saturated  poroelastic  material  is  assumed,  as  shown  in
Figure~\ref{num_model}(a). The top surface of the soil layer, $z=0$, is drained,
while its base, $z=h$, is rigid and impervious.
%
\begin{figure}[!htb]
\begin{LARGE}
\begin{sffamily}
\begin{center} 
\includegraphics[width=0.40\textwidth]{FIGS/cons_mod1.pdf}
\vspace*{-0.3cm}
\caption{\label{num_model}{(a) The geometry of the soil layer of thickness $h$ under uniform constant vertical pressure $\varpi$ = 400~kPa applied at the surface, and (b) the finite element mesh with the appropriate boundary conditions and loading at top.}}
\end{center}
\end{sffamily}
\end{LARGE}
\end{figure}


%
Since  this  is a one-dimensional problem, the only non-zero displacement is the
vertical one, $u_z$. In particular, the fluid pore pressure, $p$, as well as the displacement, $u_z$, depend only
on $z$ and $t$. Thus, the only non-zero strain component is $\varepsilon_{zz}= {\partial u_z}/{\partial z}$. 
The hydraulic boundary conditions require that $p= 0$ for the soil surface ($z= 0$) and ${\partial p}/{\partial z}=0$ for the base of the layer ($z=h$). 
The undeformability of the substratum requires that the vertical soil displacement, $u_z$, is zero ($u_z=0$ for $z=h$).
%\begin{center}
%\begin{eqnarray}
%\boldsymbol u=u(z,t)e_z \qquad p=p(z,t)
%\label{99}
%\end{eqnarray}
%\end{center}
%\begin{center}
%\begin{eqnarray}
%\varepsilon_{zz}=\frac{\partial u}{\partial z}
%\label{1099}
%\end{eqnarray}
%\end{center}
%\begin{center}
%\begin{eqnarray}
%z=0  \qquad p= 0 ; \qquad z=h  \qquad  \frac{\partial p}{\partial z}=0
%\label{95}
%\end{eqnarray}
%\end{center}
%\begin{center}
%\begin{eqnarray}
%z=h \qquad u=0
%\label{97}
%\end{eqnarray}
%\end{center}
At    time  $t=0$,  the  instantaneous  vertical  constant load $\varpi$ is suddenly applied
on the top surface $z=0$. The equilibrium requires that, for $t>0$ and $z=0$, $\sigma_{zz}=-\varpi$.
%\begin{center}
%\begin{eqnarray}
%t>0 \rightarrow \qquad z=0  \qquad \sigma_{zz}=-\varpi
%\label{96}
%\end{eqnarray}
%\end{center}
Thus, the new fields of fluid pressure and
displacement, induced by the external loading, remain to be determined.

%
Application of the constitutive Eq.~(\ref{shear_strain}) gives:
\begin{eqnarray}
\frac{\partial u_z}{\partial z}=\varepsilon_{zz}=\frac{\alpha p-\varpi}{K+(4/3)\mu}
\label{1101}
\end{eqnarray}
The instantaneous response of the porous medium to any external loading is undrained and can be expressed in the form:
\begin{eqnarray}
\varepsilon(t=0_+)=-\frac{1}{\alpha M}p(t=0_+)
\label{1103}
\end{eqnarray}
%
which, when combined with Eq.~(\ref{shear_strain}), leads to:

\begin{eqnarray}
p(z,t=0_+)
=
\frac{\alpha M\varpi}{K_u+(4/3)\mu} 
\qquad 
\frac{\partial \ud_z}{\partial z}(z,t=0_+)
=
\varepsilon_{zz}
=
\frac{-\varpi}{K_u+(4/3)\mu}
\label{1104}
\end{eqnarray}
Finally, in case of one-dimensional consolidation, the diffusion Eq.~(\ref{diffusion}), needed to be solved, reads:
\begin{eqnarray}
t>0 \qquad \frac{\partial p}{\partial t}=c_f\frac{\partial^2p}{\partial z^2}
\label{102}
\end{eqnarray}
Satisfying the above mentioned boundary conditions, the fluid pressure can be expressed in the form of infinite series:
\begin{eqnarray}
p(z,t)=\frac{\alpha M\varpi}{K_u+(4/3)\mu}\sum_{m=0}^\infty\frac{4}{\pi(2m+1)}\sin\left(\frac{(2m+1)\pi}{2}\frac{z}{h}\right)\exp\left(-\frac{(2m+1)^2\pi^2}{4}\frac{t}{\tau}\right)
\label{103}
\end{eqnarray}
Each  term  of  the  series  decreases  exponentially  with respect to the ratio
${t}/{\tau}$, with $\tau$ being the characteristic consolidation time, defined as:
\begin{eqnarray}
\tau=\frac{h^2}{c_f}
\label{104}
\end{eqnarray}
By  substituting  Eq.~(\ref{103})  into  Eq.~(\ref{1101}),  the series
converges and, after integration term by term, we obtain:
\begin{eqnarray}
u_z(z,t)  =  
\nonumber \\
\frac{\varpi(h-z)}{K+(4/3)\mu}  +  
\left(\frac{\varpi h}{K_u+(4/3)\mu}-\frac{\varpi h}{K+(4/3)\mu}\right)
\nonumber \\
\sum_{m=0}^\infty\frac{8}{\pi^2(2m+1)^2}
\cos\left(\frac{(2m+1)\pi}{2}\frac{z}{h}\right)
\exp\left(-\frac{(2m+1)^2\pi^2}{4}\frac{t}{\tau}\right)
%\nonumber \\
\label{106}
\end{eqnarray}
%
Using  Eq.~(\ref{106})  and  substitute $z=0$, the total settlement of the
soil layer, $s$, can be expressed as:
%
\begin{eqnarray}
s(t)=s_\infty+(s_{0^+}-s_\infty)\sum_{m=0}^\infty\frac{8}{\pi^2(2m+1)^2}exp\{-\frac{(2m+1)^2\pi^2}{4}\frac{t}{\tau}\} 
\end{eqnarray}
%
where
%
\begin{eqnarray}
s_{0^+}=\frac{h\varpi}{K_u+(4/3)\mu} 
\;\;\; \mbox{;} \;\;\; 
s_\infty=\frac{h\varpi}{K+(4/3)\mu}
\label{107}
\end{eqnarray}


%

{\bf Numerical Analysis} --
%
A    soil  column  of  $u$--$p$--$U$, 8-node brick, finite  elements  is used to
model     the     horizontal     layer.     Each    element    has    dimensions
1~m$\times$1~m$\times$1~m,  resulting  in  soil--column  mesh,  10  m  high,  as
illustrated in Figure~\ref{num_model}(b). The elastic material properties, shown
in  Table~\ref{Soil  Properties  for  quasi--static  problems},  are  chosen  as
indicative  values  for  the  natural  soil deposit. A constant uniform vertical
pressure  of 400~kPa is instantaneously applied on the top of the finite element
mesh.
%
\begin{table}[!htb]
\begin{center}
\caption{ \label{Soil Properties for quasi--static problems}
Soil properties used in numerical analysis of quasi--static problems.}
%\begin{spacing}{1.4}
%\setlength{\extrarowheight}{3pt}
\begin{tabular}{l|c|c|c}
       Parameter                    &  Symbol [Units]       & 1D consolidation    & Line injection                         \\
\hline
      gravity acceleration         & $g$ [m/s$^2$]        & $9.81$              & $9.81$   \\
      soil Young's Modulus         & $E$ [kN/m$^2$]       & $10^4$              & $1.2\times 10^6$ \\
      soil Poisson's ratio         & $v$                  & $0.25$              & $0.2$                        \\
      solid density                & $\rho_s$ [ton/m$^3$] & $2.65$              & $2.7$     \\
      fluid density                & $\rho_f$ [ton/m$^3$] & $1.0$               & $2.7$    \\
      solid bulk modulus           & $K_s$ [kN/m$^2$]     & $37\times 10^6$     & $3.6\times 10^7$             \\
      fluid bulk modulus           & $K_f$ [kN/m$^2$]     & $2.2\times 10^6$    & $1.0\times 10^{17}$           \\
      porosity                     & $n$                  & 0.46                & 0.4     \\
      Biot coefficient             & $\alpha$             &  $1.0$              & $1.0$\\
      Darcy permeability           & $k_D$ [m/s]          &  $10^-3$            & $3.6\times 10^{-6}$    \\
%\hline
      soil drained bulk modulus    & $K$ [kN/m$^2$]      & $6.67\times 10^3$    & $6.7\times 10^5$ \\
   soil undrained bulk modulus     & $K_u$ [kN/m$^2$]    & $4\times 10^6$       & $6.0\times 10^7$  \\
%   soil shear modulus            & $\mu$       & $4\times 10^3$~kN/m$^2$  \\
   fluid diffusivity coefficient   & $c_f$ [m$^2$/s]     & $1.2$                & $0.5$    \\
%  characteristic consolidation time    & $\tau$       & $83.3$~s  \\
\end{tabular}
%\end{spacing}
\end{center}
\end{table}




%

The    following    boundary  conditions  are  applied  to  the  model, shown in
Figure~\ref{num_model}(b):  the  solid  and fluid displacements are fixed at the
base and the pore pressure is kept zero at the top surface. In order to simulate
the  1D consolidation problem, the lateral movement of the solid and fluid phase
is  suppressed  so  that  the  vertical  displacement  remains the only non-zero
displacement.  To  remedy  any  artificial  oscillation  due to spatial and time
discretization,  numerical  damping  is  introduced  into  the analysis by using
$\gamma=0.6$    and    $\beta=0.3025$    for   the   Newmark   time   integrator
\citep{Newmark1959,local-87}.


The  instantaneous  response  of  the  soil skeleton is undrained which leads to
sudden  increase of pore pressure. The boundary condition of perfect drainage at
the  surface  ($p=0$)  causes  an  upward  flow  of the pore water, so that pore
pressure        gradually       reduces    with    time,    as   depicted   by  
Figure~\ref{dissipation}(a).  The  time factor, $T_v$ is defined as $T_v=t/\tau$
so that
for  $T_v=1$, the process of consolidation is considered to have been completed,
since      the      rate   of  dissipation  has  been  practically  diminished. 
Figure~\ref{dissipation}(b) depicts the distribution of pore pressure with depth
for
different  moments in time, $T_v$, confirming the upward movement of pore water,
as the rate of pore pressure dissipation is clearly greater for depths closer to
the  surface. 
%
\begin{figure}[!htb]
\begin{Large}
\begin{sffamily}
\begin{center}
\includegraphics[width=0.9\textwidth]{FIGS/cons_pp.pdf}
\vspace*{-0.5cm}
\caption{\label{dissipation}{Comparison    between    numerical   analysis   and
analytical  solution  by  \citet{Coussy2004}:  a)  normalized  pore  pressure at
various depths with respect to real time, $t$, and b) distribution of normalized
pore pressure with normalized depth, $z/h$, for different time factors, $T_v$.}}
\end{center}
\end{sffamily}
\end{Large}
\end{figure}
%
The gradual upward movement of pore water leads to increase of the
vertical soil displacement with time, as shown in Figure~\ref{settle}(a). It can
be concluded that the numerical analysis can effectively demonstrate the process
of  the  dissipation  of  the  pore  pressure  in  agreement with the analytical
solution by \citet{Coussy2004}.
%
\begin{figure}[tb]
\begin{Large}
\begin{sffamily}
\begin{center}
\includegraphics[width=0.4\textwidth]{FIGS/cons_disp.pdf}
\vspace*{-0.5cm}
\caption{\label{settle}{(a) Comparison between numerical analysis and analytical
solution  by  \citet{Coussy2004}  in terms of time histories: a) solid--skeleton
displacements, and b) fluid displacements, for various depths.}}
\end{center}
\end{sffamily}
\end{Large}
\end{figure}
%
Moreover,  the  numerical  analysis  based  on  $u$--$p$--$U$  formulation permits the
calculation    of    upward    fluid   movement,  $U_z$,  plotted  with  time in
Figure~\ref{settle}(b).  The upward displacement of the fluid, $U_z$, at the end
of
consolidation  has been computed as  0.39~m, more than the final downward
displacement  of  the  soil skeleton (0.33~m). Thus, an accurate estimation of
the  change  in  porosity  due  to  settlement of the soil and consequent upward
drainage  of  amount  of  pore  fluid can be performed. The amount of fluid that
escaped  is  calculated  as  $\Delta  V_f=nU_z=0.46\times  0.39\text{~m}\times  1\text{~m}\times
1\text{~m}=0.1794\text{~m}^3$.   The  change  of  solid-skeleton  volume  is  equal  to  $\Delta
V_s=(1-n)s=0.54\times  0.33\text{~m}\times  1\text{~m}\times 1\text{~m}=0.1782\text{~m}^3$. If the total initial
volume of soil is $V=10\text{~m}^3$ and the initial volume of pore fluid is $nV=4.6\text{~m}^3$,
the    final    porosity    at   the   end   of   consolidation   is   estimated
$n'={(4.6-0.1794)}/{(10-0.1782)}=0.45$,  less  than  the  initial one, $n=0.46$, indicating slight compaction of the soil layer.

%It should be mentioned
%that theory does not account for change of porosity throughout the consolidation
%process since it does not involve the fluid displacements.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\subsubsection{Line Injection of Fluid in a Reservoir}
\label{Line Injection of a fluid in a Reservoir}

%{\sc \Large I revised everything, for some reason everything was quite out of place, even eqns. Keep in mind when using the thesis you have taken this from...}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
{\bf Analytical Solution by \citet{Coussy2004}} -- 
%
A reservoir  of infinite extent composed of an isotropic, homogeneous
and  saturated  poroelastic  material is considered.  Through a cylindrical well of negligible
dimensions,  fluid mass injection  is performed in all directions
orthogonal to the well axis of symmetry, forming the $z$ axis of coordinates.
The  flow  rate of fluid mass injection is constant and equal to $q$. This  is  a  problem  of  cylindrical  symmetry, requiring the use of cylindrical
coordinates ($r$,$\theta$,$z$).  The  vectors of relative flow of fluid
mass, $\boldsymbol w$,  and solid--skeleton displacement, $\boldsymbol \xi$ are expressed as $\boldsymbol w=w(r,t)\boldsymbol {e_r}$ and $\boldsymbol \xi=\xi(r,t)\boldsymbol {e_r}$,
%
%
%\begin{center}
%\begin{eqnarray}
%\boldsymbol w=w(r,t)e_r \qquad \boldsymbol{\xi}=\xi(r)e_r
%\label{108}
%\end{eqnarray}
%\end{center}
%
where  $\boldsymbol{e_r}$  is  the unit vector along the radius. Solid--skeleton strain components are derived in the form:
\begin{eqnarray}
\varepsilon_{rr}=\frac{\partial\xi}{\partial r} 
\;\;\; \mbox{;} \;\;\; 
\varepsilon_{\theta\theta}=\frac{\xi}{r} 
\;\;\; \mbox{and} \;\;\; 
\varepsilon_{ij}=0
\label{73}
\end{eqnarray}
Thus, Eqs.~(\ref{strain_p}) and (\ref{diffusion}) are reformed in cylindrical coordinates:
\begin{eqnarray}  
\varepsilon_v=\frac{1}{r} \frac{\partial (r\xi)}{\partial r}= \frac{\alpha p}{K + (4/3)\mu}
\label{strain_p1}                
\end{eqnarray}
\begin{eqnarray}  
 \frac{\partial p}{\partial t}=c_f \frac{1}{r} \frac{\partial}{\partial r} (r \frac{\partial p}{\partial r})
\label{diffusion1}                
\end{eqnarray}
Requiring that the fluid flow reduces to zero infinitely far from the
well, $rw \rightarrow 0$ as $r \rightarrow \infty$ , the fluid mass balance relationship is applied:
\begin{eqnarray}
%rw \rightarrow 0 \qquad r \rightarrow \infty \qquad t<\infty \nonumber\\
\int_{0}^{\infty}\dot{w}_{,r}(r,t)2\pi rdr=q \qquad 
\label{110}
\end{eqnarray}
%
The solution of the set of Eqs.~(\ref{strain_p1}), (\ref{diffusion1}), (\ref{110}), and (\ref{GEconservationFlow}) for the instantaneous injection of a finite volume of fluid, $\Omega$, is given in the following form:
%
\begin{eqnarray}
p=Q\frac{K + (4/3)\mu}{K+u + (4/3)\mu}\frac{\Omega}{4\pi c_{f} t}\exp(-\frac{r^2}{4c_ft})   \nonumber
\end{eqnarray}
\begin{eqnarray}
\xi=\frac{\alpha Q\Omega}{2\pi r(K_u +(4/3)\mu)}\left[1-\exp(-\frac{r^2}{4c_ft})\right]
\label{111}
\end{eqnarray}
%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
{\bf Numerical Analysis} -- 
%
As  a  result of axisymmetry of the geometry, the model has been
developed  for  $\pi/2$  segment,  5~cm  thick,  with  radius  equal  to  1~m. A
cylindrical well is drilled at the center of the model, with 1~cm radius, having
negligible  size  when  compared  to  the  model  radius,  but  still  providing
well--shaped  elements.  
%The  final  mesh is shown in Figure~\ref{fig32}.
%%
%\begin{figure}[!htb]
%\begin{center}
%\includegraphics[width=6cm]{FIGS/Inject1.pdf}
%\vspace*{-0.5cm}
%\caption{\label{fig32} {Mesh generated for line injection problem.}}
%\end{center}
%\end{figure}
%%       
Due to
plane--strain  conditions,  all  the  movements  for  solid and fluid phases are
suppressed along the vertical axis of symmetry, $z$. The $u$ and $U$ nodes along
the  outside perimeter are fixed. It should be noted that $\Omega$ is the volume
of  the  fluid  injected  per unit of vertical well length, m$^3$/m. In order to
generate  the  volume of 1~cm$^3$/m, the corresponding fluid displacement of the
nodes  along  the well has been calculated and applied as a step function at the
time  $t=0$~sec.  For  simplicity, the initial fluid pressure $p_0$ is set to be
0~kPa.  The  analytical  solution  is  studied below using the set of parameters
shown in Table~\ref{Soil Properties for quasi--static problems}.



In  the analysis, the pore pressure and the radial displacements are recorded at
three  points  on  the  radius  (10~cm,  50~cm,  and  85~cm). The analytical and
numerical      results    compare    well  with  each  other  as  indicated  by 
Figure~\ref{fig33}. 
%
\begin{figure}[~htb]
\begin{center}
\includegraphics[width=6cm]{FIGS/Injection.pdf}
%\includegraphics[width=8cm]{FIGS/InjectionPressure_new_ps1.pdf}
%\hfill
%\includegraphics[width=7.8cm]{FIGS/InjectionDisplacement_new_ps1.pdf}
\vspace*{-0.3cm}
\caption{\label{fig33}   {Comparison between numerical analysis and analytical solution by \citet{Coussy2004}: a) radial displacement, and b) pore pressure versus time at different distances, R, from the center of the model.}}
\end{center}
\end{figure}
%
The time step, $\Delta t$, was set equal to 1~sec; thus, the
first data point is recorded at the time of 1~sec.




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\clearpage
\subsection{Poroelastodynamics}  
\label{Poroelastodynamics}


One  of  the  significant    findings   by \citet{Biot1956a,Biot1956b,Biot1956c} was the
identification of two kinds of compressional waves and one shear wave, which has
been  the  starting  point  for  further research on the mechanics of dispersive
waves  through  analytical  or  numerical  approaches.  Several  analytical  and
semi-analytical  solutions  for  dynamic  poroelasticity  are  available  in the
literature,  varying  with geometry, type of loading, time domain or transformed
domain  solution,  assumptions  related  to the compressibility of the fluid and
the formulation used ($u$--$p$, $u$--$p$--$U$, $u$--$U$).
%
A  review  of  this  wide range of solutions in dynamic poroelasticity is nicely
presented   by  \citet{Schanz2009}.  The  current  paper  focuses  on  solutions
regarding  one-dimensional  wave  propagation in saturated linear elastic porous
media due to step and sinusoidal excitations.


In  order  to study the accuracy of $u$--$p$--$U$ finite element procedure in
representing   a  dispersive  wave  propagation  in  saturated  soil,  numerical
dispersion effects should be taken into account. Dependence of the wave-velocity
on  the  frequency  content  of  the  excitation is the basic concept behind the
theory  of  dispersive  waves.  Thus,  when  wave  propagation  is  modeled in a
numerical  scheme,  numerical  features,  such  as the element size, integration
scheme, numerical damping, element type etc., should be chosen appropriately, so
that  the  frequency  range  of the excitation will be represented in the model.
Otherwise,  the  wave-velocities  developed  in the numerical model will diverge
from  the  theoretical  solution  due  to an artificial dispersion attributed to
numerics and not physics (numerical dispersion).

Various   theoretical   works   concern   the  effect  of  numerical  dispersion
\citep{Deraemaeker1999,local-86,Ihlenburg1995,Semblat2000},       highlighting
spatial  discretization  as  one  of the most influential numerical features. In
particular,  the  effect of numerical dispersion decreases with the increase in wave length,
$\lambda$,  and  increases  with the increase in mesh size $\Delta h$, which can lead to the
necessity  to use fine meshes, even uneconomical ones for small wave lengths
\citep{Deraemaeker1999}. A common rule of thumb for typical finite element
procedures  is  to  resolve the wavelength by 10 elements \citep{Ihlenburg1995, local-86};
thus  the  element  size,  $\Delta  h$, should be less that 1/10 of the smallest
wavelength ($\lambda_{\min}$)  involved in the application:

\begin{eqnarray}
\Delta h \le  \frac{\lambda_{\min}}{10} 
\mbox{~~~~~~~~~where~~~~~~~~} 
\lambda_{\min}= \frac{V_{\min}}{f_{\max}}
\label{110a}
\end{eqnarray}
%
where  $V_{\min}$  is  the  lowest  wave  velocity  of  interest in the
simulation and generally, it is a function of the elastic modulus of the medium,
while  $f_{\max}$  is  the  maximum  frequency that is modeled.


Spatial  and  temporal  discretization  are  linked in a way requiring that mesh
spacing  is  accompanied  by  an  appropriate  time  interval.  As  a wave front
progresses  in  space  it reaches one point after the other. If the time step in
the  finite  element  analysis  is  too  large,  the  wave  front  can appear to reach two
consecutive  nodes  at  the  same  moment.  This  would violate a fundamental
property  of  wave  propagation  and can lead to instability. In order to ensure
that the propagation of  wave can be modeled properly, the required time
step, $\Delta t$, must be limited to:
%
\begin{eqnarray}
\Delta t \le  \frac{\Delta h}{V_{\max}}
\label{111a}
\end{eqnarray}
%
where  $V_{\max}$ is the highest wave velocity involved, usually being a function
of the undrained elastic modulus in case of fully saturated porous media.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\subsubsection{Step displacement excitation}
%\label{One dimensional shock wave propagation with a step displacement boundary condition at the surface}
\label{1DDisp_BC}
%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
{\bf Analytical Solution by \citet{Gajo1995b}} -- 
%
An  one-dimensional  analytical solution of the Biot's equations (using $u$--$U$
field  displacements)  is  provided by \citet{Gajo1995b} for a general transient
problem  in  poroelasticity.  The analytical solution was obtained using Fourier
series  expansion  and  it  is  not based on any assumptions with respect to the
inertial,  viscous  or mechanical coupling. Furthermore, it is applicable to any
type of boundary-initial value problem.

Since  each  term  of the Fourier series represents a frequency component of the
excitation  signal,  the  analytical solution can describe the behavior of every
single  frequency  component related to the problem. Thus, it can illustrate the
mechanics  of  dispersive wave propagation in which higher frequencies propagate
with  two  waves  and  lower  frequencies  with  only one wave, as a function of
permeability  and  travel length. Considering the above, the analytical solution
can  provide  a useful comparative tool towards the verification of the existing
numerical solutions based on the finite element method \citep{Gajo1994}.

%

\citet{Gajo1995b} present the response of the porous medium only relative to the
first  arrival  of  the  waves  of  the  first  and  second  kind. Specifically,
analytical  results  are  shown  for  a  finite  soil column subjected to a step
displacement boundary condition (Heaviside function) at the surface applied both
to solid and fluid phases, as shown in Figure~\ref{Fig1001}(a). 
%
\begin{figure}[!htb]
\begin{Large}
\begin{sffamily}
\begin{center}
\includegraphics[width=0.40\textwidth]{FIGS/Gajo_mod.pdf}
\caption{\label{Fig1001}{a)  A  representative  soil  column subjected to a step
vertical  displacement equal to $1.0\times10^{-3}\text{~cm}$ at the surface, and
b)  the  finite  element  mesh  and  the  applied  boundary  conditions used for
numerical modeling.}}
\end{center}
\end{sffamily}
\end{Large}
\end{figure}
%
This problem can
demonstrate  better the mechanics of dispersive wave propagation, since the step
excitation  contains waves of all frequencies. In the following section (3.2.2),
analytical  solution  is  compared  with  numerical  analysis  for a soil column
subjected  to  a  single sine pulse with characteristic predominant frequency of
$10^6\text{~Hz}$, producing smoother wavefronts.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


{\bf Numerical Analysis} -- 
%
Fine  spatial  discretization  is  required due to step boundary condition. This
kind  of  excitation  contains  waves of all frequencies leading to  meshes with
a  large  number  of finite elements, according  to Eq.~(\ref{110a}), which
suggests that the appropriate
element  size  is  inversely proportional to the maximum frequency introduced to
the  model.  In  fact, there is no element size that could satisfy the equation;
thus,  numerical  dispersion  is  inevitable.  It  can  only  be  diminished  by
increasing  the  spatial  refinement.  Evidently,  a compromise needs to be made
between the quantity of numerical error and the computational cost.

Numerical  models  with  three  different values of permeability -- resulting in
three   different   viscous  coupling  scenarios  ($k=10^{-8}\text{~cm$^3$s/g}$,
$k=10^{-6}\text{~cm$^3$s/g}$,  $k = 10^{-5}\text{~cm$^3$s/g}$) -- were solved in
order to verify the $u$--$p$--$U$ formulation in a wide range of drag forces. At
the   top   surface   of   the   model,   a   vertical   step   displacement  of
$1.0\times10^{-3}\text{~cm}$  is applied both to the solid and the fluid phases.
No  lateral  flow or displacement is allowed. The base of the model is rigid and
impervious, as schematically illustrated in Figure~\ref{Fig1001}(a).

According  to  \citet{Gajo1995},  the  two  extreme values of wave velocity for
longitudinal  propagation  in a porous medium with properties shown in the third
column  of  Table~\ref{Soil Properties_gajo}, are: $V_{\min}=610\text{~m/s}$ and
$V_{\max}=1819\text{~m/s}$,  corresponding  to  the  case  of  significantly low
viscous  coupling,  where  two longitudinal waves exist (one slow and one fast).
%

\begin{table}[!htb]
\begin{center}
\caption{ \label{Soil Properties_gajo}
Soil Properties for 1D wave propagation due
to various boundary conditions and loading.}
%\begin{spacing}{1.4}
\begin{tabular}{l|c|c|c|c}
       Parameter                   &  Symbol [Units]       &  Step/Sine pulse disp.     &  Step/Sine loading     &  Step vel.        \\
\hline
      gravity acceleration         & $g$ [m/s$^2$]         & $9.81$                            & $9.81$                       & $9.81$                \\
      soil Young's Modulus         & $E$ [kN/m$^2$]        & $1.2\times 10^6$                  & $20\times 10^3$              & $23.2\times 10^6$     \\
      soil Poisson's ratio         & $v$                   & $0.3$                             & $0.2$                        & $0.17$                \\
      solid density                & $\rho_s$ [ton/m$^3$]  & $2.7$                             & $2.0$                        & $2.66$                \\
      water density                & $\rho_f$ [ton/m$^3$]  & $1.0$                             & $1.0$                        & $1.0$                 \\
      solid bulk modulus           & $K_s$ [kN/m$^2$]      & $36\times 10^6$                   & $36\times 10^6$              & $36\times 10^6$       \\
      fluid bulk modulus           & $K_f$ [kN/m$^2$]      & $2.2\times 10^6$                  & $2.2\times 10^6$             & $2.2\times 10^6$      \\
      porosity                     & $n$                   & $0.4$                             & $0.33$                       & $0.18$                \\
      Biot coefficient             & $\alpha$              & $1.0$                             & $1.0$                        & $0.677$               \\
\end{tabular}
%\end{spacing}
\end{center}
\end{table}
%
For  infinitely  large viscous coupling, the velocity of the single longitudinal
wave  is  $1773\text{~m/s}$.  The  chosen  finite  element  mesh consists of 400
$u$--$p$--$U$   brick  finite  elements  of  dimensions  $0.01\text{~cm}  \times
0.01\text{~cm}  \times  0.01\text{~cm}$,  creating  a  soil column 4 cm high, as
illustrated  in  Figure~\ref{Fig1001}(b).  The  maximum  frequency allowed to be
present   in   a   model   with  the  above-mentioned  grid  spacing  is  $18.19
\times{10^6}\text{~Hz}$   for  the  first  longitudinal  wave  (fast)  and  $6.1
\times{10^6}\text{~Hz}$ for the second one (slow). There is a cut-off of greater
frequencies imposed by the specific computational mesh.




The    maximum    time    step    is    equal   to
$(0.0001\text{~m})/(1819\text{~m/s})=5.5\times{10^{-8}}\text{~sec}$,          as
suggested  by  Eq.~(\ref{111a}).  The selected temporal integration involves 800
steps  of  $\Delta t = 2.0\times10^{-8}\text{~sec}$, which allows a maximum wave
velocity   of   $5.0\times10^5\text{~m/s}$.  
% In   addition,  the  Newmark  time
%integration  method  was  used \citep{Newmark1959}.
% Sets of parameters, assuming
%unconditionally numerical stability, were chosen for both integrators.


Figure~\ref{Fig13}  shows the comparative results for all three different values
of  viscous  coupling  using  Newmark  parameters:  $\gamma  = 0.6$ and $\beta =
0.3025$,  which introduce  numerical damping in the model. 
%
\begin{figure}[!htb]
\begin{Large}
\begin{sffamily}
\begin{center}
\includegraphics[width=0.7\textwidth]{FIGS/Gajo_comp.pdf}
\vspace*{-0.5cm}
\caption{\label{Fig13}{  Comparison  between  numerical  analysis and analytical
solution  by  \citet{Gajo1995b}  for  three different values of viscous coupling,
$k$:  a)  solid  displacement versus  time  from 4--16$\times$10$^{-6}$~sec, b) solid
displacement  versus  time  from  5--6.5$\times$10$^{-6}$~sec, c) fluid displacement versus
time  from  4--16$\times$10$^{-6}$~sec, and  d)  fluid  displacement versus time from
5--6.5$\times$10$^{-6}$~sec, obtained $1$~cm below the surface.}}
\end{center}
\end{sffamily}
\end{Large}
\end{figure}
%  
In  general,  the  numerical  results    are   in good agreement with analytical
solution  for  1D wave propagation in fully saturated, elastic porous media. For
example,  numerical  analysis  demonstrates  that  during the propagation of the
first  wave,  the  solid  and  fluid  displacement are in phase with each other,
whereas during the propagation of second wave, the displacements are in opposite
phase.  In  case  of  high  viscous coupling, only one longitudinal wave exists,
since the relative motion between solid and fluid phase is highly constrained.

However,  as expected,  the rise time of wave fronts obtained from the numerical
analysis  is  larger than the rise time estimated by the analytical solution, due
to  numerical  dispersion,  observed also by  \citet{Gajo1994}  and
\citet{Simon1986}.  This  discrepancy is smaller for high viscous coupling which
physically  induces  a  more  dispersive  response,  in contrast to low
viscous  coupling which produces  abrupt  wavefronts.  This  dispersive  response in case of large viscous coupling is
attributed  to  the  physically  based, high viscous damping due to coupling of pore
fluid  and  porous solid, which can be captured by the $u$--$p$--$U$ formulation
(see matrices $C_1$, $C_2$, and $C_2$ of Eq.~(\ref{69-1}), in the appendix).

Furthermore, the rise time of the second wave is even larger; a fact that can be
partly attributed to the diffusive behavior of the second wave in low--frequency
range.  In  general,  the  distortion  and  the smearing of the wave fronts is a
drawback  of all types of numerical solutions, which is linked to the limitation
of  the  highest  frequency  allowed by the computational grid and the numerical
damping introduced to the model.


In  order  to  further  examine the sensitivity of the numerical solution to
characteristic numerical  features,  such  as  a)  spatial refinement and b) numerical
damping, an indicative parametric investigation was conducted.  For  this  purpose,  a  finer  spatial  discretization  of  $\Delta h=
0.005\text{~cm}$  was  selected,  accompanied  by  a  smaller time step equal to
$10^{-8}\text{~sec}$.  Two  different  sets  of parameters for Newmark integrator
were  used:  i)  $\gamma = 0.5$ and $\beta = 0.25$ corresponding to no numerical
damping   and  ii) $\gamma = 0.6$ and $\beta = 0.3025$, which introduce numerical
damping. Figure~\ref{NM05} shows
the  numerical  response versus the analytical solution for two combinations
of  spatial  and  temporal refinement when no numerical damping is introduced in
the  model  ($\gamma  =  0.5$  and  $\beta  =  0.25$). 
%
\begin{figure}[!htb]
\begin{Large}
\begin{sffamily}
\begin{center}
\includegraphics[width=0.7\textwidth]{FIGS/gajo_step_NM05.pdf}
\vspace*{-0.3cm}
\caption{\label{NM05}{ Analytical solution by \citet{Gajo1995b} versus numerical
analysis  with  no numerical damping ($\gamma = 0.5$ and $\beta = 0.25$) for two
different   combinations  of  spatial  and  temporal  discretization:i)  $\Delta
h=0.01\text{~cm}$    and    $\Delta t = 2 \times 10^{-8}\text{~sec}$,    ii)  
$\Delta
h=0.005\text{~cm}$  and  $\Delta  t=10^{-8}\text{~sec}$,  in  the case of viscous
coupling $k = 10^{-5}\text{~cm$^3$s/g}$.}}
\end{center}
\end{sffamily}
\end{Large}
\end{figure}
%
High-frequency artificial
oscillations  develop  within the model at the first wave arrival with gradually
decreasing  amplitude  over  time. 
%
%
%
Evidently,  the oscillations are narrower with equal  or  smaller  amplitude  in
case  of  the  finer  mesh, while the rise time of abrupt wavefronts diminishes,
better approaching the analytical solution.

Figure~\ref{NM06}  shows  the  comparison  of  numerical results with analytical
solution  for  the  two  selected  meshes  when numerical damping is used in the
numerical  scheme  ($\gamma  =  0.6$  and  $\beta  = 0.3025$). 
%
\begin{figure}[!htb]
\begin{Large}
\begin{sffamily}
\begin{center}
\includegraphics[width=0.7\textwidth]{FIGS/gajo_step_NM06.pdf}
\vspace*{-0.3cm}
\caption{\label{NM06}{ Analytical solution by \citet{Gajo1995b} versus numerical
analysis with no numerical damping ($\gamma = 0.6$ and $\beta = 0.3025$) for two
different   combinations  of  spatial  and  temporal  discretization:i)  $\Delta
h=0.01\text{~cm}$    and    $\Delta    t=2x10^{-8}\text{~sec}$,    ii)    $\Delta
h=0.005\text{~cm}$  and  $\Delta  t=10^{-8}\text{~sec}$,  in  the case of viscous
coupling $k = 10^{-5}\text{~cm$^3$s/g}$.}}
\end{center}
\end{sffamily}
\end{Large}
\end{figure}
%
Once more, closer
convergence  with  analytical  solution is achieved in case of finer mesh. It is
interesting  to  notice  that  the  slight  oscillations  occurring  close to the
wavefronts  in  case  of  coarser  mesh,  have  been completely damped out by
numerical  damping  in  case of finer spatial discretization. This fact leads to
the   conclusion  that  numerical  damping  is, in this case,   more  effective 
in  damping  out  higher-frequency  oscillations  which develop in case of finer
mesh.

Overall,  Figure~\ref{zoom}  suggests  that  the  analytical  solution is better
approached in terms of the rise time of wavefronts, when no numerical damping is
introduced  in  the  model  and finer mesh is used. 
%
\begin{figure}[!htb]
\begin{Large}
\begin{sffamily}
\begin{center}
\includegraphics[width=0.7\textwidth]{FIGS/gajo_step_zoom.pdf}
\vspace*{-0.3cm}
\caption{\label{zoom}{ Analytical solution by \citet{Gajo1995b} versus numerical
analysis  for  two  cases  of  numerical damping: a) $\gamma = 0.5$ and $\beta =
0.25$  (no numerical damping) and b) $\gamma = 0.6$ and $\beta = 0.3025$ and for
two  different  combinations  of spatial and temporal discretization: i) $\Delta
h=0.01\text{~cm}$  and  $\Delta  t=2  \times  10^{-8}\text{~sec}$,  ii)  $\Delta
h=0.005\text{~cm}$   and  $\Delta  t=10^{-8}\text{~sec}$,  in  case  of  viscous
coupling  $k = 10^{-5}\text{~cm$^3$s/g}$. Magnified view of the first wavefront.
}}
\end{center}
\end{sffamily}
\end{Large}
\end{figure}
%
However,  numerical  damping  is  needed     to    damp    out  the  artificial 
oscillations  developed  due  to  the  discretization  of  continuum into finite
element  model  in  combination  with  the  range  of frequencies imposed by the
excitation.
%
The    compromise    in   choosing the parameters of the numerical features  is 
directly  related  to  the desired level of accuracy in capturing the details of
response, where both low and high frequency waves are present.
% It should be noted that high
%frequency waves  are damped out faster.

%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\subsubsection{Single sine pulse displacement}
\label{1--D  wave propagation due to single sine pulse displacement}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

In order to investigate the efficiency of $u$--$p$--$U$ finite element under more
realistic excitations with a characteristic predominant frequency, a single sine
pulse  displacement  was  imposed to both solid and fluid phases at the top of a
finite  soil  column.  The  characteristic  duration of the single sine pulse is
$10^{-6}\text{~sec}$,    corresponding    to    a   predominant   frequency   of
$10^6\text{~Hz}$.  The  soil properties assigned to the model are the same as in
the  previous  example  and  are  shown in the third column of Table~\ref{Soil Properties_gajo}. The
minimum   wavelength,   $\lambda_{\min}$,   is   estimated   as  $V_{\min}/f_{\max}=
(610\text{~m/s})/(10^6\text{Hz})   =   6.1   \times  10^{-4}\text{~m}$.  According  to
Eq.~(\ref{110a}),  the  element  size is  limited to $1/10$ of the minimum
wavelength,   so   that   it   is  less  than  $6.1  \times{10^{-5}}\text{~m}$. 
The selected mesh is comprised of $800$ $u$--$p$--$U$ brick finite elements     
of  dimensions  $0.005\text{~cm}\times{0.005\text{~cm}}\times{0.005\text{~cm}}$,
creating  a  soil column $4\text{~cm}$ high. The time step is limited to 
$\Delta t \le \Delta
h/V_{\max}  =  2.7 \times{10^{-8}}\text{~sec}$ according to Eq.~(\ref{111a}). The
chosen  temporal  discretization  consists  of time intervals equal to 
$\Delta t = 2 \times 10^{-8}\text{~sec}$.   
%
Numerical   models   with   three   different  values  of
permeability, resulting  in  three  different  viscous  coupling cases ($k =
10^{-8}\text{~cm$^3$s/g}$,   $k  =  5  \times  10^{-8}\text{~cm$^3$s/g}$,  and  $k  =
10^{-5}\text{~cm$^3$s/g}$), were analyzed.




Figure~\ref{gajo_sine}  shows  the  comparative  results for all three different
values of viscous coupling using Newmark parameters: $\gamma = 0.5$ and $\beta =
0.25$,  introducing  no  numerical  damping.  
%
\begin{figure}[!htb]
\begin{Large}
\begin{sffamily}
\begin{center}
\includegraphics[width=0.7\textwidth]{FIGS/gajo_sine.pdf}
\vspace*{-0.3cm}
\caption{\label{gajo_sine}{  Comparison  between  numerical  analysis and analytical
solution    by    \citet{Gajo1995b}  in  terms  of  solid displacement, obtained
$1\text{~m}$  below the surface, for three different values of viscous coupling,
$k$.}}
\end{center}
\end{sffamily}
\end{Large}
\end{figure}
%
%
%
In  general,  the  numerical  results  converge  well  with analytical solution,
especially  for  higher  values  of  viscous  coupling  where  the  response  is
physically  more  dispersive.  One  can  observe  that  the  there  is  a slight
discrepancy  between analytical and numerical response in terms of the rise time
at  the  arrival  of  the  pulse, for cases of lower viscous coupling. This fact emphasizes the
sensitivity  of  the  response  to  the  spatial  refinement even under smoother
excitations  with a predominant characteristic frequency, particularly for cases
of  low viscous coupling. The numerical error can be eliminated with finer mesh,
as  shown  in  the  previous  case.  However,  herein,  the  effect of numerical
dispersion  is significantly smaller compared to the example of step excitation and it could be suppressed in case of completely smooth pulse (e.g. Ricker).
Moreover,  no  need  for  numerical  damping  occurred  since  no high frequency
artificial oscillations developed as in case of step pulse.


%old
%In  general, the numerical results
%compare  well with analytical solution, especially for higher values of viscous
%coupling  where the response is physically more dispersive. One can observe that
%the  there  is discrepancy between analytical and numerical response in 
%the rise time at the arrival of the pulse due to numerical dispersion. This fact
%emphasizes  the sensitivity of the response to the spatial refinement even under
%smoother  excitations,  particularly  for  cases  with low viscous coupling. The
%numerical  error  can  be  eliminated  with finer mesh, as shown in the previous
%case.  There  was  no  need for numerical damping since there were no artificial
%high frequency oscillations as in case of step loading.
%



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{Step loading}
\label{One dimensional shock wave propagation with step loading at the surface}

{\bf Analytical Solution by \citet{deBoer1993}} --
%
An      analytical     solution    for    one-dimensional  wave  propagation  in
fluid-saturated elastic porous media is provided by \citet{deBoer1993}, based on
theory  of  mixtures \citep{Bowen1980,Bowen1982,Ehlers1993}. The fluid-saturated
porous  material  is assumed to be a two phase system composed of incompressible
solid  and  fluid  phases.  An exact analytical solution is obtained via Laplace
transform technique exhibiting only one independent compressive wave in both the
solid and fluid phases, as a result of the incompressibility constraint.

The  problem  configuration consists  of  an infinitely  long column,  extracted
from  the  half-space of a fluid-saturated porous elastic skeleton material. The
motions  of  both solid and fluid materials are constrained to occur only in the
vertical direction. Loading is applied to the surface boundary, as a function of
time,  $\sigma(z=0,t)  =  f(t)$.  Both solid and fluid displacements, as well as
solid skeleton extra stress and pore pressure are obtained with respect to time and
depth \citep{deBoer1993}. The analytical solution has been
developed    for    i) sinusoidal, ii) step  (presented here) and iii) impulsive
loading.

{\bf Numerical Analysis} -- 
%
Numerical  example  for  the  step  loading  case  was  used to verify the
implementation  of  $u$--$p$--$U$  formulation. Step loading excitation contains
all frequencies, thus requiring a fine spatial discretization. However, this
example  does  not  examine  the details of wavefront; it rather focuses on the long
term  response  describing  the  consolidation  process  under  dynamic loading.
Considering  the  above,  the  numerical  grid used for the simulation of the 1D
shock  wave  propagation consists of 1000 $u$--$p$--$U$ brick finite elements of
dimensions   $1\text{~cm}\times1\text{~cm}\times1\text{~cm}$   creating  a  soil
column  $10$~m  high.  Since  the  numerical  simulation of a semi infinite soil
column  is  not possible, a soil column of finite depth of $10$~m was considered
adequate  for the current problem configuration. No lateral flow or displacement
are allowed in the model; thus, only the vertical displacement is free. The pore
pressure  is  constrained  at  the  top  surface  to be equal to the atmospheric
pressure; thus, modeling drained condition. The base of the model is rigid and impervious.

Table~\ref{Soil Properties_gajo}  shows  the  soil properties of the numerical 
model,  which  are the same with those used for the analytical results presented
in  the  paper by \citet{deBoer1993}. The only difference is associated with the
elastic  modulus,  which  was  selected  to  be  $20$~MN/m$^2$ for the numerical
solution instead of $30$~MN/m$^2$, as noted by \citet{deBoer1993}. The value of Darcy permeability assigned to the soil skeleton is equal to 0.01 m/s. The numerical
analysis indicates that the analytical results presented in the paper correspond
to  a  soil  column  with  elastic  modulus  equal  to  $20$~MN/m$^2$ instead of
$30$~MN/m$^2$. Moreover, it should be mentioned that the solid and fluid (water)
compressibility    were    given    realistic    values   (see  the fourth column of Table~\ref{Soil Properties_gajo}),  which  practically  means  that  the two constituents can be
assumed fairly incompressible compared to the soil--skeleton compressibility.


At    the  top  surface  of  the soil column, a step loading of $\sigma(z=0,t) =
3\text{~kN/m$^2$}$  is applied to the solid part through four nodal loads at the
top. The nodal loads are developed as:
%
\begin{eqnarray}
 F_N(z=0,t) 
 = 
 \frac{\sigma(z=0,t)\times{A}}{4} 
 = 
 \frac{3\text{~kN/m$^2$}\times0.01\text{~m}\times0.01\text{~m}}{4} 
 = 
 7.5\times10^{-5}\text{~kN}
\label{1}
\end{eqnarray}

The maximum wave velocity of the soil can be given approximately by equation:
%
\begin{eqnarray}
V_{\max}
=
\sqrt{\frac{M+\alpha\times{Q}}{(1-n)\times{\rho_s}+n\times{\rho_f}}}
=
2000\text{~m/s} 
\label{vmax}
\end{eqnarray}
%
where    $M$  is the constrained modulus of the solid skeleton, $\alpha$ is Biot
coefficient, while $Q$ is the bulk stiffness of a mixture, defined after Eq.
(\ref{GEconservationFlow}) on page \pageref{GEconservationFlow}.
%    %
Thus,   according  to  Eq.~(\ref{111a})  the  time  step  should  be 
$5\times{10^{-6}}\text{~sec}$ or less.  The  selected temporal integration involves 80000
steps  of $\Delta t = 5.0\times10^{-6}\text{~sec}$. The Newmark parameters used in the analysis are: $\gamma = 0.7$ and $\beta =
0.4$.  This  sets of parameters introduce considerable numerical damping to the
model  which  can  increase  numerical  dispersion, but as mentioned before, the
current  study  does not focuses on the details of the high frequency components
which are quickly damped out with time.
%

Figures~\ref{Fig_disp},  \ref{Fig_stress},  and \ref{Fig_pp}  show the
comparison  between  analytical and numerical results. 
%
\begin{figure}[!htb]
\begin{Large}
\begin{sffamily}
\begin{center}
\includegraphics[width=0.7\textwidth]{FIGS/Boer_disp.pdf}
\vspace*{-0.5cm}
\caption{\label{Fig_disp}{Comparison   between numerical analysis and analytical
solution by \citet{deBoer1993} in terms of a) solid displacements versus time at
different  depths,  b)  fluid  displacements versus time at different depths, c)
solid  displacements versus depth at different times, and d) fluid displacements
versus depth at different times.}}
\end{center}
\end{sffamily}
\end{Large}
\vspace{-.5cm}
\end{figure}

\begin{figure}[!htb]
\begin{Large}
\begin{sffamily}
\begin{center}
\includegraphics[width=0.7\textwidth]{FIGS/Boer_stress.pdf}
\vspace*{-0.5cm}
\caption{\label{Fig_stress}{Comparison      between   numerical   analysis   and
analytical  solution by \citet{deBoer1993} in terms of a) solid--skeleton stress
versus  time  at different depths, and b) solid--skeleton stress versus depth at
different times.}}
\end{center}
\end{sffamily}
\end{Large}
\vspace{-.5cm}
\end{figure}

\begin{figure}[!htb]
\begin{Large}
\begin{sffamily}
\begin{center}
\includegraphics[width=0.7\textwidth]{FIGS/Boer_pp.pdf}
\vspace*{-0.5cm}
\caption{\label{Fig_pp}{Comparison    between  numerical analysis and analytical
solution  by \citet{deBoer1993}: a) pore fluid pressure versus time at different
depths, b) pore fluid pressure stress versus depth at time $t=0.01$~sec, c) pore
fluid  pressure  stress  versus  depth  at  time $t=0.05$~sec, and d) pore fluid
pressure stress versus depth at time $t=0.1$~sec.}}
\end{center}
\end{sffamily}
\end{Large}
\vspace{-.5cm}
\end{figure}
%
In general, the numerical
results  are  in  good  agreement,  with  respect  to time and depth, with those
obtained by the analytical solution. During the consolidation process, the solid
skeleton  settles and the fluid escapes out of the solid skeleton while the load
is  transferred  from  the  pore pressure to the solid skeleton stress. The only
discrepancy  between  analytical  and  numerical  results  is  found in the pore
pressures  which develop an oscillatory trend in case of numerical analysis. The
existence  of  oscillatory  waves  has  also  been  observed  and  discussed  by
\citet{Zienkiewicz84}.   These   artificial   oscillations,  attributed  to  the
combination  of  numerical  features and the nature of dynamic step loading, are
damped  out  over time. Numerical damping would be more efficient in eliminating
the  artificial  oscillations  if finer mesh was used and thus, higher-frequency
oscillations  were developed.  This  conclusion,  indicating  that  the  efficiency of
numerical  damping  increases as the artificial oscillations become narrower, was
shown in section \ref{1DDisp_BC} and will be also shown in section \ref{1DVelBC}.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\subsubsection{Sinusoidal loading}
\label{1--D wave propagation due to sinusoidal loading}


Numerical  example  for the sinusoidal loading case was also used to verify 
implementation  of $u$--$p$--$U$ formulation. The frequency of the excitation is
equal    to   $75$   rad/s   and   the   load   is   given   by   the   function
$3\times{(1-\text{cos}(75\times{t}))}$.  Table~\ref{Soil Properties_gajo} shows
the  soil  properties  of  the  numerical  model,  which  are the same as in the
previous section. The value of Darcy permeability assigned to the soil skeleton is again equal to 0.01 m/s.


The  minimum wave velocity that can be present in the model can be approximately
estimated by equation:
%
\begin{eqnarray}
V_{\min}=\sqrt{\frac{M}{(1-n)\times{\rho_s}}}=128\text{~m/s} 
\label{vmin}
\end{eqnarray}
%
corresponding  to  a  dry  solid  skeleton.  Thus, the minimum wavelength can be
estimated  by  Eq.~(\ref{110a})  as $10.7\text{~m}$ leading to a maximum element
size  of $1.07\text{~m}$. The numerical grid used for the simulation of the 1--D
shock  wave  propagation  consists of 100 $u$--$p$--$U$ brick finite elements of
dimensions  $10\text{~cm} \times 10\text{~cm} \times 10\text{~cm}$  creating  a
soil  column  $10\text{~m}$  high. The maximum time-step required is estimated by
Eq.~(\ref{111a})  equal  to  $4\times{10^{-6}}\text{~sec}$ approximately for the
selected  element size. The Newmark parameters were selected to be: $\gamma  =  0.6$  and $\beta = 0.3025$, in order to introduce numerical damping. At the top
surface  of the soil column, the sinusoidal loading is applied to the solid part
through four nodal loads at the top.

Figures~\ref{Boer_sine}   depicts   the   time-histories   of  solid  and  fluid
displacements, as well as pore pressure and effective stress as obtained by both
numerical  analysis and analytical solution. 
%
\begin{figure}[!htb]
\begin{Large}
\begin{sffamily}
\begin{center}
\includegraphics[width=0.7\textwidth]{FIGS/deBoer_sine.pdf}
\vspace*{-0.5cm}
\caption{\label{Boer_sine}{Comparison  between numerical analysis and analytical
solution  by  \citet{deBoer1993}  in  terms  of a) solid displacements, b) fluid
displacements,  c)  pore  pressure  and d) solid skeleton stress versus time for
different depths.}}
\end{center}
\end{sffamily}
\end{Large}
\vspace{-.5cm}
\end{figure}
%
Numerical analysis captures
the response of fully saturated soil under sinusoidal
loading   in a satisfactory manner.   In  particularly,  the  solid  displacement  gradually
accumulates downwards
while  the  fluid is squeezed out of the skeleton. However, there are moments of
partial recovery, when the load on the surface decreases, which renders the pore
pressure  values  negative  close  to  surface, as the fluid is absorbed inwards
(suction).  The  effective  stress  seems to be more sensitive to loading close to the
surface  following  the  shape  of  the  excitation, whereas the effect tends to
diminish for greater depths due to solid-fluid interaction.



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


\subsubsection{Step velocity excitation}
%\label{One dimensional shock wave propagation with a step velocity boundary condition}
\label{1DVelBC}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
{\bf Semi--Analytical Solution by \citet{Hiremath88}} -- 
%
\citet{Hiremath88}   present  a semi--analytical solution of Biot's equations of
motion  for one-dimensional wave propagation in a fluid-saturated linear elastic
isotropic  medium  (soil)  using  Laplace  transformation  followed by numerical
inversion.  This  study  is considered to be an extension of the exact transient
solution  presented by \citet{Garg74} for two limiting cases of infinitely small
and  infinitely  large  viscous coupling. In both cases, a soil column of finite
dimension  (50~cm)  subjected to velocity boundary conditions at the surface was
analyzed, allowing for reflection of waves at the boundaries.

One    of   the   most  important  observations,  made by  both  \citet{Garg74} 
and  \citet{Hiremath88}, is that in case of strong viscous coupling (high drag),
the material behaves as a single continuum with internal dissipation and the two
wave fronts tend to become a single one.

In    particular,   \citet{Hiremath88}  examined  two  cases  allowing for small
and  large  viscous  coupling. Moreover, two different types of excitations were
applied  at  the  boundary  surface in terms of solid and fluid velocity. In the
first  case,  a  unit step boundary condition was applied at the top surface for
both solid and fluid phases (to be studied below). In the second case, the fluid
velocity    applied    at    the   boundary  is slightly different  from  the
applied solid
velocity,  and  it is increasing gradually to unity over a short time scale. The
results obtained
from   the   numerical  inversion  allowed  for  six  reflections  of  the  fast
compressional  wave  of  first  kind  and  two reflections of the secondary slow
longitudinal wave.



{\bf Numerical Analysis} -- 
%
Models      for      two     extreme    values  of  viscous  coupling:  a)  large
($k=0.148\times10^{-8}\text{~cm$^3$s/g}$)           and          b)          low
($k=0.148\times10^{-2}\text{~cm$^3$s/g}$)  are  developed and simulated. At the
top surface of the soil column, a step velocity of $1.0\times10^{-2}\text{~m/s}$
is   applied  both  to  the  solid  and  the  fluid  phase.  Only  the  vertical
translational   degrees  of  freedom  are  free  so  that  no  lateral  flow  or
displacement   are  allowed.  The  base of the model is rigid and impervious. A 
step excitation introduces all frequencies in the model,
requiring  a very small  element  size (theoretically, based on Eq.~(\ref{110a}),
the  element  size  should really be infinitely small, in order to accommodate all
introduced frequencies). The
selected finite element mesh consists of 100 $u$--$p$--$U$ brick finite elements
of  dimensions $0.005\text{~m}\times0.005\text{~m}\times0.005\text{~m}$ creating
a  soil  column  50~cm high, as suggested by \citet{Hiremath88}. Table~\ref{Soil Properties_gajo}  shows  the  soil properties used for this model. The solid
skeleton  is so rigid that the wave of the first kind propagates with a velocity
close  to 3500m/s according to equation Eq.~(\ref{vmax}). The required time step
is     estimated     by     Eq.~(\ref{111a})   and is  equal    to    $\Delta    t    =
1.4\times10^{-6}\text{~sec}$.  The  selected  temporal integration involves 1972
steps of $\Delta t = 5.0\times10^{-7}\text{~sec}$. The Newmark set of parameters
was  selected  as  $\gamma  =  0.6$  and $\beta = 0.3025$, introducing numerical
damping in the model.



Figures~\ref{solid_10_high}    and    \ref{solid_30_high}    show the comparison
between  numerical  and  analytical  results  for  both extreme cases of viscous
coupling.  
%
\begin{figure}[!htb]
\begin{Large}
\begin{sffamily}
\begin{center}
\includegraphics[width=0.7\textwidth]{FIGS/Hiremath_high.pdf}
\vspace*{-0.5cm}
\caption{
\label{solid_10_high} 
Comparison   between   numerical   analysis   and   analytical   solution   
\citet{Hiremath88}     for     the    case    of    strong    viscous    coupling
($k=0.148\times10^{-8}\text{~cm$^3$s/g}$)    in  terms of time histories of: a) solid velocity at
$10$~cm  below  the  surface,  b) solid velocity at $30$~cm below the surface, c) fluid
velocity  at  $10$~cm  below the surface and, d) fluid velocity at $30$cm below the
surface.}
\end{center}
\end{sffamily}
\end{Large}
\end{figure}

\begin{figure}[!htb]
\begin{Large}
\begin{sffamily}
\begin{center}
\includegraphics[width=0.7\textwidth]{FIGS/Hiremath_low.pdf}
\vspace*{-0.5cm}
\caption{\label{solid_30_high}{Comparison   between   numerical   analysis   and
analytical  solution  \citet{Hiremath88}  for  the  case of weak viscous coupling
($k=0.148\times10^{-2}\text{~cm$^3$s/g}$):  a) solid
velocity  at  $10$~cm  below  the  surface,  b) solid velocity at $30$~cm below the
surface,  c)  fluid velocity at $10$~cm below the surface and d) fluid velocity at
$30$cm below the surface, versus time.}}
\end{center}
\end{sffamily}
\end{Large}
\end{figure}
%
Overall,  the finite element solution reproduces correctly the trends
of  wave  propagation in both limiting cases of viscous coupling. In particular,
numerical  analysis  demonstrates  well  that  for  the  case  of strong viscous
coupling,  solid and fluid velocities are in phase with each other. 
%It should be
%noted  that  the  numerical  analysis develops oscillations at the instantaneous
%changes in velocity due to reflection of wave fronts.


In  order to improve the accuracy of the numerical solution, a finer spatial and
temporal discretization was used, with $\Delta h = 0.001\text{~m}$ and $\Delta t
=  2\times{10^{-7}}\text{~sec}$,  while  Newmark  parameters  remained  the same.
%
Figure~\ref{hiremath_NM06}  depicts the analytical solution versus the numerical
results  obtained from both meshes. 
%
\begin{figure}[!htb]
\begin{Large}
\begin{sffamily}
\begin{center}
\includegraphics[width=0.7\textwidth]{FIGS/hiremath_NM06.pdf}
\vspace*{-0.5cm}
\caption{\label{hiremath_NM06}{Analytical   solution  \citet{Hiremath88}  versus
numerical  analysis,  with  two  different  combinations of spatial and temporal
discretization  and numerical damping ($\gamma = 0.6$ and $\beta = 0.3025$), for
the  case of weak viscous coupling ($k=0.148\times10^{-2}\text{~cm$^3$s/g}$): a)
solid  velocity  at 10 cm below the surface and b) fluid velocity at 10 cm below
the surface, with time.}}
\end{center}
\end{sffamily}
\end{Large}
\end{figure}
% 
Evidently, the finer mesh achieves a greater
accuracy,  suppresses oscillations,      and  gives a more similar response  to the  analytical 
solution, as expected.
%
It  is interesting to notice that when numerical damping is used, all artificial
oscillations are suppressed in case of finer mesh, while they persist in case of
coarser mesh, as illustrated in Figure~\ref{hiremath_fine}.
%

\begin{figure}[!htb]
\begin{Large}
\begin{sffamily}
\begin{center}
\includegraphics[width=0.7\textwidth]{FIGS/hiremath_fine.pdf}
\caption{\label{hiremath_fine}{Analytical   solution  \citet{Hiremath88}  versus
numerical  analysis  with  no  numerical  damping  ($\gamma  = 0.6$ and $\beta =
0.3025$),  in  terms  of fluid velocity time history at 10 cm below the surface,
for the case of weak viscous coupling ($k=0.148\times10^{-2}\text{~cm$^3$s/g}$):
a) $\Delta h = 0.001\text{~m}$ and $\Delta t = 2 \times 10^{-7}\text{~sec}$ and 
b) $\Delta h = 0.005\text{~m}$ and $\Delta t = 5 \times 10^{-7}\text{~sec}$.}}
\end{center}
\end{sffamily}
\end{Large}
\end{figure}

       


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%\clearpage

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Conclusions}

A  solution verification suite for modeling and
simulation  of fully coupled behavior of saturated porous media, using u-p-U formulation,
have been presented in this paper. 
%
Analytical solutions for both static and
dynamic  examples have been included in this suite. 
%
%
Detailed  comparison  between numerical  and analytical solutions shows   close 
matching,  while  common  numerical  discretization  effects  are  captured  and
controlled.   In   particular,   artificially--introduced   (by  discretization)
higher-frequency  effects  have been observed, as in most numerical simulations.
These  numerical  dispersion  effects  have  been  investigated  focusing on the
influence of numerical features such as spatial and temporal refinement, as well
as  numerical damping on the obtained response. It is illustrated that the finer
meshes  better  approach  the  analytical  solutions, as expected; however, such
finer meshes produce even higher frequency (artificial) oscillations. It is also
shown  that  such  oscillations produced by finer meshes can be more efficiently
damped out numerically than those produced by coarser discretization.


%
Presented solution verification suite can be used by any 
modeling and simulation effort that deals with the response of fully saturated porous 
media. In this particular case, presented solution verification results provide 
evidence that the fully coupled models for
saturated porous media, from the UCD Computational Geomechanics Group numerical 
libraries, are solved correctly for both static and 
dynamic cases.


It  is  important  to  observe  that  any  numerical method develop to model and
simulation  behavior  of  solids/structures and fluids needs to be verified (and
validated   if   quality  experimental  data  is  available).  Only  after  such
verification  (and validation) process has been successfully completed, a method
can  be  claimed  to  be  complete.  The  main aim of this paper, with presented
verification  of  the  u-p-U  formulation,  is to contribute to such approach to
numerical modeling and simulations.




%Presented solution verification results  provide evidence that
%the fully coupled model for saturated porous media was solved correctly by UCD  Computational  Geomechanics   
%Group's
%implementation    of  3D  $u$--$p$--$U$    finite  elements \citep{Jeremic2008}. The developed solution verification procedure can be used for any finite element or finite difference code oriented to simulate the response of fully saturated porous media.

%The next step would be 
%developing  a  validation  suite, providing evidence that the correct model (or set of
%models) is (are) solved. For validation purposes, high quality experimental results are
%required; thus, finding the most appropriate ones, is our next goal.


%\clearpage
\begin{appendix}

%\section{APPENDIX}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Finite Element Formulation}   
\label{Finite Element Formulation}


Standard finite element discretization \citep{FEM-IV-1,FEM-IV-2}, using shape
functions, is used to describe 
each  of  the  unknown  fields  ($u$--$p$--$U$) in terms of nodal values 
(solid  displacements  $\bar{u}_{Ki}$, pore fluid pressures $\bar{p}_K$, and pore
fluid displacements $\bar{U}_{Ki}$):
%
\begin{equation}                 
  u_i = N_K^u \bar{u}_{Ki}, \;\;\;             
  p = N_K^p \bar{p}_K, \;\;\;            
  U_i = N_K^U \bar{U}_{Ki}           
\label{39}                 
\end{equation}    
%            
where    $N_K^u$,   $N_K^p$, and $N_K^U$ are (in this case same) shape functions
for solid displacement,
pore    pressure,  and  fluid  displacement,  respectively.  Each  node  of  the
($u$--$p$--$U$)  element  thus features seven degrees of freedoms in three dimensions
(three  for solid displacements, one for pore fluid pressures, and three or pore
fluid  displacements).  It should be noted that it is possible to use same shape
functions  for  both  displacement  and  pore  pressure  unknown  field  as  the
$u$--$p$--$U$ formulation with compressible fluid allows that without volumetric
locking.


By using finite element discretization and  after   some   tensor  algebra  and 
manipulations, the weak form of governing equations can be
obtained  from  the  strong  form  described by Eqs.~(\ref{34}), (\ref{35}), and
(\ref{36}):
%                                                                                                                                             %%
%%
%
%
%
%
\begin{eqnarray}
& &                                                                                                                                         
\left[ \begin{array}{ccc}  
(M_s)_{KijL} & 0 & 0 \\                                                                                 
 0  & 0 & 0 \\  
 0  & 0 & (M_f)_{KijL}
\end{array} \right]
\left[ \begin{array}{c}
\ddot{\bar{u}}_{Lj} \\
\ddot{\bar{p}}_N \\
\ddot{\bar{U}}_{Lj}
\end{array} \right]    
+   
\left[ \begin{array}{ccc}
(C_1)_{KijL} & 0 & -(C_2)_{KijL} \\
 0  & 0 &   0  \\  
-(C_2)_{LjiK} & 0 & (C_3)_{KijL} \\     
\end{array} \right]
\left[ \begin{array}{c}
\dot{\bar{u}}_{Lj} \\
\dot{\bar{p}}_N \\
\dot{\bar{U}}_{Lj}
\end{array} \right] \nonumber \\      
&+&    
\left[ \begin{array}{ccc}   
0 & -(G_1)_{KiM} & 0 \\                                                                                                                            
-(G_1)_{LjM}  & -P_{MN} & -(G_2)_{LjM} \\  
 0  & -(G_2)_{KiL} & 0
\end{array} \right]
\left[ \begin{array}{c}   
\bar{u}_{Lj} \\
\bar{p}_M \\
\bar{U}_{Lj}
\end{array} \right] \nonumber \\
&+& 
\left[ \begin{array}{c}   
\int_{\Omega} N_{K,j}^u \sigma^{'}_{ij} \ud V \\                                                      
0  \\  
0  
\end{array} \right]     
-   
\left[ \begin{array}{c}   
\bar{f}^u_{Ki} \\                                                      
0  \\  
\bar{f}^U_{Ki}  
\end{array} \right] 
= 0    
\label{68}      
\end{eqnarray}
%%     
where the left hand side components are:                                                                                                                                            %%
% %%%  
%%%  
\begin{eqnarray}
(M_s)_{KijL}= \int_{V} N_K^u (1-n) \rho_s \delta_{ij} N_L^u \ud V
\;\;\;\; &\mbox{;}& \;\;\;\;
(M_f)_{KijL}= \int_{V} N_K^U n \rho_f \delta_{ij} N_L^U \ud V  \nonumber\\  
(C_1)_{KijL}= \int_{V} N_K^u n^2 k_{ij}^{-1} N_L^u \ud v 
\;\;\;\;   &\mbox{;}& \;\;\;\;
(C_2)_{KijL}= \int_{V} N_K^u n^2 k_{ij}^{-1} N_L^U \ud V  \nonumber\\  
(C_3)_{KijL}= \int_{V} N_K^U n^2 k_{ij}^{-1} N_L^U \ud V 
\;\;\;\; &\mbox{;}& \;\;\;\;
(G_1)_{KiM} = \int_{V} N_{K,i}^u (\alpha-n) N_M^p \ud V   \nonumber\\  
(G_2)_{KiM} = \int_{V} n N_{K,i}^U N_M^p \ud V    
\;\;\;\; &\mbox{;}& \;\;\;\;
P_{NM}      = \int_{V} N_N^p \frac{1}{Q} N_M^p \ud V 
\label{69-1}  
\end{eqnarray}
%%%
%
%
while the right hand side are:
% %%%   
\begin{eqnarray}                                                                       
(\bar{f}_s)_{Ki} &=& \int_{A_t} N_K^u \sigma^{'}_{ij} n_j \ud A  
                     - 
                     \int_{A_p} N_K^u (\alpha-n) p n_i \ud A  
                     + 
                     \int_{V} N_K^u (1-n) \rho_s b_i \ud V  \nonumber\\                     
(\bar{f}_f)_{Ki} &=& - \int_{A_p} N_K^U n p n_i \ud A  
                     + 
                     \int_{V} N_K^U n \rho_f b_i \ud V   
\label{69-2}  
\end{eqnarray}  
%
where  $V$  is  the  volume, and $A_t$ and $A_p$ are the domain boundaries with
traction and the pore fluid pressure defined, respectively.
 
It        is    very   important  to note that the velocity proportional viscous
damping  was  introduced  directly  through  the  damping tensor with components
$(C_1)_{KijL}$,  $(C_2)_{KijL}$  and  $(C_3)_{KijL}$,  which  are  functions  of
porosity  and permeability of the skeleton. This damping provides for physically
based, velocity proportional energy dissipation due to interaction of pore fluid
and  the  solid  (soil)  skeleton.  It  is also emphasized 
that presented
formulation and implementation do not (need to) use Rayleigh damping.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 \section{Time Integration}   
 \label{Time_Integration}
 
 
 In  order  to develop integration of dynamic finite element equation in the time
 domain,  Eq.~(\ref{68})  is rewritten in a residual matrix form \citep{local-87}:
 %%%
 \begin{eqnarray}   
   \bf{ R = M \ddot{x} + C \dot{x} + K' x + F(x) - f = 0 }          
 \label{71} 
 \end{eqnarray}
 %%%
 where  $\bf{x  =  \{\bar{u},  \bar{p},  \bar{U}\}^T}$  represent  a  vector  of 
 generalized unknown variables.
 Eq.~(\ref{71}) represents the general non-linear form  \citep{local-87} for which 
 the usual tangent stiffness $\bf K$ obtained from:  
 \begin{eqnarray}    
   \bf{ K = \frac{\partial R}{\partial x} = K' + \frac{\partial F(x)}{\partial x} }                                
 \label{72}    
 \end{eqnarray}
 
 In  the  specific  case  of  $u$--$p$--$U$  formulation of interest here one can write
 matrix form:
 %%%
 \begin{eqnarray}    
   \bf{ M  = 
     \left[ \begin{array}{ccc}       
        \bf{M_s} & 0 & 0 \\               
        0  & 0 & 0 \\        
        0  & 0 & \bf{M_f}             
     \end{array} \right] }, \;\;\;                 
   \bf{ C  = 
      \left[ \begin{array}{ccc}       
        \bf{C_1} & 0 & \bf{-C_2} \\                  
        0  & 0 &   0  \\      
       \bf{-C_2^T} & 0 & \bf{C_3} \\                   
      \end{array} \right]  }                                                                                                               
 \label{upU_MC}        
 \end{eqnarray}
 %%
 \begin{eqnarray}
  \bf{ K'  =                                                                                                                                                
   \left[ \begin{array}{ccc}       
      0 & \bf{-G_1} & 0 \\       
      \bf{-G_1^T}  & \bf{-P} & \bf{-G_2^T}\\                     
      0  & \bf{-G_2} & 0           
    \end{array} \right]},\;\;\;              
   \bf{ K  = 
   \left[ \begin{array}{ccc}      
      \bf{K^{ep}} & \bf{-G_1} & 0 \\                 
      \bf{-G_1^T}  & \bf{-P} & \bf{-G_2^T}\\                   
      0  & \bf{-G_2} & 0         
    \end{array} \right] }                                                                                                                
 \label{upU_K}      
 \end{eqnarray}
 where
 \begin{eqnarray}
   {\bf K^{ep}} = (K^{ep})_{KijL} = \int_{V} N_{K,m}^u E^{ep}_{imjn} N_{L,n}^u \ud V                                                
 \label{upU_Kep}      
 \end{eqnarray}
 
 
 The  above  set  of  residual  (nonlinear)  dynamic  equations  is solved using the  
 Newmark   procedure \citep{Newmark1959}. This time integration method has two parameters,
 $\beta$ and $\gamma$, and is described by the following two equations:
 %%%%
 %{\sc \Large I checked this, same with Lecture Notes}
 \begin{eqnarray}
   {}^{n+1} x &=& {}^n x + \Delta t \; {}^n \dot{x} + \Delta^2 t \; [(\frac{1}{2} - \beta) \; {}^n \ddot{x} + \beta \; {}^{n+1} \ddot{x}]
 \label{upU_Newmark_EQ1} \\
   {}^{n+1} \dot{x} &=& {}^n \dot{x} + \Delta t \; [(1 - \gamma) \; {}^n \ddot{x} + \gamma \; {}^{n+1} \ddot{x}] 
 \label{upU_Newmark_EQ2}      
 \end{eqnarray}
 %%%%
 which give the relations between the current time step $n$ to the next time step $n+1$. 
 The method is generally an implicit method, except when both $\beta$ and $\gamma$ are zero. 
%
%{\sc \Large I checked this, same with Lecture Notes}
%
If the parameters $\beta$ and $\gamma$ satisfy the following conditions
%
 \begin{equation}
   \gamma \ge \frac{1}{2}, \;\;\; \beta = \frac{1}{4}(\gamma+\frac{1}{2})^2
 \label{upU_Newmark_ab}     
 \end{equation}
 %
 the  time  integration  method  is  unconditionally  stable.  Any $\gamma$ value
 greater  than  $0.5$ will introduce numerical damping. Well-known members of the
 Newmark  time  integration  method  family  include: trapezoidal rule or average
 acceleration  method  for  $\beta = 1/4$ and $\gamma = 1/2$, linear acceleration
 method  for  $\beta=1/6$  and  $\gamma = 1/2$, and (explicit) central difference
 method  for  $\beta  = 0$ and $\gamma = 1/2$. If and only if  $\gamma = 1/2$, it is second order accurate \citep{local-86}.
% 
% 
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% 
% \section{Program Implementation}   
% \label{PI}
% 
% %{\sc \Large NEEDS TO BE REWRITTEN}
% 
% Implementation  of  the  described  algorithms  and  procedures  was performed using a
% number    of    numerical    libraries,  algorithms,  material models, and finite
% elements  available through one of open source licenses.
% %
% %
% The finite  element  domain  is  managed  using  Modified OpenSees Services (MOSS)
% library 
% \citep{Jeremic2004d,  Jeremic2007d,  Jeremic2008a,   McKenna97}.
% %
% On      a   lower functional  level,   a   set   of   Template3Dep   numerical   libraries
% \citep{Jeremic2000f}
% are used  for  constitutive  level modeling, while nDarray numerical libraries
% \citep{Jeremic97d} are used to handle vector, matrix
% and  tensor  manipulations, and FEMtools  element  libraries from UCD
% CompGeoMech      toolset    \citep{Jeremic2004d}    are  used to supply other
% necessary libraries and components.
% %
% %
% %
% The solution    of  system of equations is provided by  UMFPACK solver
% \citep{Davis97b, Davis97, Davis1999, Davis2004b, Davis2004}.
% %




\end{appendix}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%\include{Appendix}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%$
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%$


%\clearpage
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\bibliography{refmech,Mahdi,zhao}
%\bibliography{./refmech,./refcomp,./Mahdi,./zhao,./ref_Giota}
\bibliography{refmech,refcomp,Mahdi,zhao,ref_Giota}
\bibliographystyle{abbrvnat}
%\bibliographystyle{agsm}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



       
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%$
%TABLES
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%$
% \clearpage
% \begin{table}[!htb]
% \begin{center}
% \caption{ \label{Soil Properties for quasi--static problems}
% Soil properties used in numerical analysis of quasi--static problems.}
% %\begin{spacing}{1.4}
% %\setlength{\extrarowheight}{3pt}
% \begin{tabular}{l|c|c|c}
%        Parameter                    &  Symbol [Units]       & 1D consolidation    & Line injection                         \\
% \hline
%       gravity acceleration         & $g$ [m/s$^2$]        & $9.81$              & $9.81$   \\
%       soil Young's Modulus         & $E$ [kN/m$^2$]       & $10^4$              & $1.2\times 10^6$ \\
%       soil Poisson's ratio         & $v$                  & $0.25$              & $0.2$                        \\
%       solid density                & $\rho_s$ [ton/m$^3$] & $2.65$              & $2.7$     \\
%       fluid density                & $\rho_f$ [ton/m$^3$] & $1.0$               & $2.7$    \\
%       solid bulk modulus           & $K_s$ [kN/m$^2$]     & $37\times 10^6$     & $3.6\times 10^7$             \\
%       fluid bulk modulus           & $K_f$ [kN/m$^2$]     & $2.2\times 10^6$    & $1.0\times 10^{17}$           \\
%       porosity                     & $n$                  & 0.46                & 0.4     \\
%       Biot coefficient             & $\alpha$             &  $1.0$              & $1.0$\\
%       Darcy permeability           & $k_D$ [m/s]          &  $10^-3$            & $3.6\times 10^{-6}$    \\
% %\hline
%       soil drained bulk modulus    & $K$ [kN/m$^2$]      & $6.67\times 10^3$    & $6.7\times 10^5$ \\
%    soil undrained bulk modulus     & $K_u$ [kN/m$^2$]    & $4\times 10^6$       & $6.0\times 10^7$  \\
% %   soil shear modulus            & $\mu$       & $4\times 10^3$~kN/m$^2$  \\
%    fluid diffusivity coefficient   & $c_f$ [m$^2$/s]     & $1.2$                & $0.5$    \\
% %  characteristic consolidation time    & $\tau$       & $83.3$~s  \\
% \end{tabular}
% %\end{spacing}
% \end{center}
% \end{table}
% 

%\begin{table}[!htb]
%\begin{center}
%\caption{\label{table8}{Soil properties used to study the line injection problem.}}
%\begin{spacing}{1.4}
%\setlength{\extrarowheight}{3pt}
%\begin{tabular}{l|c|c}
%      Parameter                    & Symbol    &  Value                         \\
%\hline
%Young's modulus & E & $1.2\times 10^6$ kN/m$^2$ \\
%Poisson ratio & $\nu$ & $0.2$  \\ 
%solid density & $\rho_s$ & $2.7$ ton/m$^3$ \\
%fluid density & $\rho_f$ & $1.0$ ton/m$^3$ \\
%solid bulk modulus & $K_s$ & $3.6\times 10^7$ kN/m$^2$ \\
%fluid bulk modulus & $K_f$ & $1.0\times 10^{17}$ kN/m$^2$ \\
%porosity & $n$ & 0.4  \\
%Darcy permeability & $k_{D}$ & $3.6\times 10^{-6}$ m/s \\
%bulk modulus & $K$ & $6.7\times 10^5$ kN/m$^2$ \\
%undrained bulk modulus & $K_u$ & $6.0\times 10^7$ kN/m$^2$ \\
%fluid diffusivity coefficient & $c_f$ & $0.4973$ m$^2$/s \\
%\end{tabular}
%\end{spacing}
%\end{center}
%\end{table}


% \begin{table}[!htb]
% \begin{center}
% \caption{ \label{Soil Properties_gajo}
% Soil Properties for 1D wave propagation due
% to various boundary conditions and loading.}
% %\begin{spacing}{1.4}
% \begin{tabular}{l|c|c|c|c}
%        Parameter                   &  Symbol [Units]       &  Step/Sine pulse displacement     &  Step/Sinusoidal loading     &  Step velocity        \\
% \hline
%       gravity acceleration         & $g$ [m/s$^2$]         & $9.81$                            & $9.81$                       & $9.81$                \\
%       soil Young's Modulus         & $E$ [kN/m$^2$]        & $1.2\times 10^6$                  & $20\times 10^3$              & $23.2\times 10^6$     \\
%       soil Poisson's ratio         & $v$                   & $0.3$                             & $0.2$                        & $0.17$                \\
%       solid density                & $\rho_s$ [ton/m$^3$]  & $2.7$                             & $2.0$                        & $2.66$                \\
%       water density                & $\rho_f$ [ton/m$^3$]  & $1.0$                             & $1.0$                        & $1.0$                 \\
%       solid bulk modulus           & $K_s$ [kN/m$^2$]      & $36\times 10^6$                   & $36\times 10^6$              & $36\times 10^6$       \\
%       fluid bulk modulus           & $K_f$ [kN/m$^2$]      & $2.2\times 10^6$                  & $2.2\times 10^6$             & $2.2\times 10^6$      \\
%       porosity                     & $n$                   & $0.4$                             & $0.33$                       & $0.18$                \\
%       Biot coefficient             & $\alpha$              & $1.0$                             & $1.0$                        & $0.677$               \\
% \end{tabular}
% %\end{spacing}
% \end{center}
% \end{table}
% 

%\begin{table}[!htb]
%\begin{center}
%\caption{ \label{Soil Properties_boer}
%Soil Properties for 1D wave propagation due
%to step loading.}
%\begin{spacing}{1.4}
%\begin{tabular}{l|c|c}
%      Parameter                    &  Symbol    &  Value                         \\
%\hline
%      gravity acceleration         & $g$       & 9.81~m/s$^2$         \\
%      soil matrix Young's Modulus  & $E$       & $20\times 10^3$~kN/m$^2$              \\
%      soil matrix Poisson's ratio  & $v$       & 0.2                              \\
%      solid particle density       & $\rho_s$  & $2.0\times 10^3$~kg/m$^3$        \\
%      water density                & $\rho_f$  & $1.0\times 10^3$~kg/m$^3$        \\
%      solid particle bulk modulus  & $K_s$     & $36.0\times 10^6$~kN/m$^2$              \\
%      fluid bulk modulus           & $K_f$     & $2.177\times 10^6$~kN/m$^2$              \\
%      porosity                     & $n$       & 0.33                             \\
%      Darcy permeability         & $k_D$     & 0.01~m/s                         \\
%\end{tabular}
%\end{spacing}
%\end{center}
%\end{table}

%\begin{table}[!htb]
%\begin{center}
%\caption{ \label{Soil Properties_Hiremath}
%Soil Properties for 1D wave propagation due to step velocity boundary condition.}
%\begin{spacing}{1.4}
%\begin{tabular}{l|c|c}
%       Parameter                    &  Symbol    &  Value                         \\
%\hline
%      gravity acceleration         & $g$       & 9.81~m/s$^2$         \\
%      soil matrix Young's Modulus  & $E$       & $23.21\times 10^6$~kN/m$^2$              \\
%      soil matrix Poisson's ratio  & $v$       & 0.171                              \\
%      solid particle density       & $\rho_s$  & $2.66\times 10^3$~kg/m$^3$        \\
%      fluid density                & $\rho_f$  & $1.0\times 10^3$~kg/m$^3$        \\
%      solid particle bulk modulus  & $K_s$     & $36.0\times 10^6$~kN/m$^2$              \\
%      fluid bulk modulus           & $K_f$     & $2.2\times 10^6$~kN/m$^2$              \\
%      porosity                     & $n$       & 0.18                            \\
%      Biot coefficient             & $\alpha$  &  $0.6772$                    \\
%\end{tabular}
%\end{spacing}
%\end{center}
%\end{table}



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%$
%FIGURES
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%$
% \clearpage
% 
% %
% \begin{figure}[!htb]
% \begin{center}
% {\includegraphics[width=9.0cm]{FIGS/RoleVV_NEW_knowledge.pdf}}
% \caption{\label{RoleVV} Schematic representation of the role of Verification and
% Validation  on  numerical  modeling  (inspired     from     \protect{\citet{Oberkampf2002,Roy2011}}     and
% \protect{\citet{Oden2010a}}).}
% \end{center}
% \end{figure}



% \begin{figure}[!htb]
% \begin{LARGE}
% \begin{sffamily}
% \begin{center} 
% \includegraphics[width=0.50\textwidth]{FIGS/cons_mod1.pdf}
% \caption{\label{num_model}{(a) The geometry of the soil layer of thickness $h$ under uniform constant vertical pressure $\varpi$ = 400~kPa applied at the surface, and (b) the finite element mesh with the appropriate boundary conditions and loading at top.}}
% \end{center}
% \end{sffamily}
% \end{LARGE}
% \end{figure}


% \begin{figure}[!htb]
% \begin{Large}
% \begin{sffamily}
% \begin{center}
% \includegraphics[width=1\textwidth]{FIGS/cons_pp.pdf}
% \caption{\label{dissipation}{Comparison    between    numerical   analysis   and
% analytical  solution  by  \citet{Coussy2004}:  a)  normalized  pore  pressure at
% various depths with respect to real time, $t$, and b) distribution of normalized
% pore pressure with normalized depth, $z/h$, for different time factors, $T_v$.}}
% \end{center}
% \end{sffamily}
% \end{Large}
% \end{figure}

% \begin{figure}[!htb]
% \begin{Large}
% \begin{sffamily}
% \begin{center}
% \includegraphics[width=0.5\textwidth]{FIGS/cons_disp.pdf}
% \caption{\label{settle}{(a) Comparison between numerical analysis and analytical
% solution  by  \citet{Coussy2004}  in terms of time histories: a) solid--skeleton
% displacements, and b) fluid displacements, for various depths.}}
% \end{center}
% \end{sffamily}
% \end{Large}
% \end{figure}

%        
% \begin{figure}[!htb]
% \begin{center}
% \includegraphics[width=8cm]{FIGS/Injection.pdf}
% %\includegraphics[width=8cm]{FIGS/InjectionPressure_new_ps1.pdf}
% %\hfill
% %\includegraphics[width=7.8cm]{FIGS/InjectionDisplacement_new_ps1.pdf}
% \caption{\label{fig33}   {Comparison between numerical analysis and analytical solution by \citet{Coussy2004}: a) radial displacement, and b) pore pressure versus time at different distances, R, from the center of the model.}}
% \end{center}
% \end{figure}
       
% \begin{figure}[!htb]
% \begin{center}
% \includegraphics[width=8cm]{FIGS/Inject1.pdf}
% \caption{\label{fig32} {Mesh generated for line injection problem}}
% \end{center}
% \end{figure}
       

% \begin{figure}[!htb]
% \begin{Large}
% \begin{sffamily}
% \begin{center}
% \includegraphics[width=0.50\textwidth]{FIGS/Gajo_mod.pdf}
% \caption{\label{Fig1001}{a) A representative soil column subjected to a step vertical displacement equal to $1.0\times10^{-3}\text{~cm}$ at the surface, and b) the finite element mesh and the applied boundary conditions used for numerical modeling.}}
% \end{center}
% \end{sffamily}
% \end{Large}
% \end{figure}
       
       
% \begin{figure}[!htb]
% \begin{Large}
% \begin{sffamily}
% \begin{center}
% \includegraphics[width=1\textwidth]{FIGS/Gajo_comp.pdf}
% \caption{\label{Fig13}{  Comparison  between  numerical  analysis and analytical
% solution  by  \citet{Gajo1995b}  for  three different values of viscous coupling,
% $k$:  a)  solid  displacement versus  time  from 4--16$\times$10$^{-6}$~sec, b) solid
% displacement  versus  time  from  5--6.5$\times$10$^{-6}$~sec, c) fluid displacement versus
% time  from  4--16$\times$10$^{-6}$~sec, and  d)  fluid  displacement versus time from
% 5--6.5$\times$10$^{-6}$~sec, obtained $1$~cm below the surface.}}
% \end{center}
% \end{sffamily}
% \end{Large}
% \end{figure}

% \begin{figure}[!htb]
% \begin{Large}
% \begin{sffamily}
% \begin{center}
% \includegraphics[width=1\textwidth]{FIGS/gajo_step_NM05.pdf}
% \caption{\label{NM05}{ Analytical solution by \citet{Gajo1995b} versus numerical
% analysis  with  no numerical damping ($\gamma = 0.5$ and $\beta = 0.25$) for two
% different   combinations  of  spatial  and  temporal  discretization:i)  $\Delta
% h=0.01\text{~cm}$    and    $\Delta t = 2 \times 10^{-8}\text{~sec}$,    ii)  
% $\Delta
% h=0.005\text{~cm}$  and  $\Delta  t=10^{-8}\text{~sec}$,  in  the case of viscous
% coupling $k = 10^{-5}\text{~cm$^3$s/g}$.}}
% \end{center}
% \end{sffamily}
% \end{Large}
% \end{figure}


% \begin{figure}[!htb]
% \begin{Large}
% \begin{sffamily}
% \begin{center}
% \includegraphics[width=1\textwidth]{FIGS/gajo_step_NM06.pdf}
% \caption{\label{NM06}{ Analytical solution by \citet{Gajo1995b} versus numerical
% analysis with no numerical damping ($\gamma = 0.6$ and $\beta = 0.3025$) for two
% different   combinations  of  spatial  and  temporal  discretization:i)  $\Delta
% h=0.01\text{~cm}$    and    $\Delta    t=2x10^{-8}\text{~sec}$,    ii)    $\Delta
% h=0.005\text{~cm}$  and  $\Delta  t=10^{-8}\text{~sec}$,  in  the case of viscous
% coupling $k = 10^{-5}\text{~cm$^3$s/g}$.}}
% \end{center}
% \end{sffamily}
% \end{Large}
% \end{figure}

% \begin{figure}[!htb]
% \begin{Large}
% \begin{sffamily}
% \begin{center}
% \includegraphics[width=1\textwidth]{FIGS/gajo_step_zoom.pdf}
% \caption{\label{zoom}{ Analytical solution by \citet{Gajo1995b} versus numerical
% analysis  for  two  cases  of  numerical damping: a) $\gamma = 0.5$ and $\beta =
% 0.25$  (no numerical damping) and b) $\gamma = 0.6$ and $\beta = 0.3025$ and for
% two  different  combinations  of spatial and temporal discretization: i) $\Delta
% h=0.01\text{~cm}$  and  $\Delta  t=2  \times  10^{-8}\text{~sec}$,  ii)  $\Delta
% h=0.005\text{~cm}$   and  $\Delta  t=10^{-8}\text{~sec}$,  in  case  of  viscous
% coupling  $k = 10^{-5}\text{~cm$^3$s/g}$. Magnified view of the first wavefront.
% }}
% \end{center}
% \end{sffamily}
% \end{Large}
% \end{figure}

% \begin{figure}[!htb]
% \begin{Large}
% \begin{sffamily}
% \begin{center}
% \includegraphics[width=1\textwidth]{FIGS/gajo_sine.pdf}
% \caption{\label{gajo_sine}{  Comparison  between  numerical  analysis and analytical
% solution    by    \citet{Gajo1995b}  in  terms  of  solid displacement, obtained
% $1\text{~m}$  below the surface, for three different values of viscous coupling,
% $k$.}}
% \end{center}
% \end{sffamily}
% \end{Large}
% \end{figure}

% \begin{figure}[!htb]
% \begin{Large}
% \begin{sffamily}
% \begin{center}
% \includegraphics[width=1\textwidth]{FIGS/Boer_disp.pdf}
% \caption{\label{Fig_disp}{Comparison  between numerical analysis and analytical solution by \citet{deBoer1993} in terms of a) solid displacements versus time at different depths, b) fluid displacements versus time at different depths, c) solid displacements versus depth at different times, and d) fluid displacements versus depth at different times.}}
% \end{center}
% \end{sffamily}
% \end{Large}
% \vspace{-.5cm}
% \end{figure}
% 
% \begin{figure}[!htb]
% \begin{Large}
% \begin{sffamily}
% \begin{center}
% \includegraphics[width=1\textwidth]{FIGS/Boer_stress.pdf}
% \caption{\label{Fig_stress}{Comparison  between numerical analysis and analytical solution by \citet{deBoer1993} in terms of a) solid--skeleton stress versus time at different depths, and b) solid--skeleton stress versus depth at different times.}}
% \end{center}
% \end{sffamily}
% \end{Large}
% \vspace{-.5cm}
% \end{figure}
% 
% \begin{figure}[!htb]
% \begin{Large}
% \begin{sffamily}
% \begin{center}
% \includegraphics[width=1\textwidth]{FIGS/Boer_pp.pdf}
% \caption{\label{Fig_pp}{Comparison  between numerical analysis and analytical solution by \citet{deBoer1993}: a) pore fluid pressure versus time at different depths, b) pore fluid pressure stress versus depth at time $t=0.01$~sec, c) pore fluid pressure stress versus depth at time $t=0.05$~sec, and d) pore fluid pressure stress versus depth at time $t=0.1$~sec.}}
% \end{center}
% \end{sffamily}
% \end{Large}
% \vspace{-.5cm}
% \end{figure}
% 

% \begin{figure}[!htb]
% \begin{Large}
% \begin{sffamily}
% \begin{center}
% \includegraphics[width=1\textwidth]{FIGS/deBoer_sine.pdf}
% \caption{\label{Boer_sine}{Comparison  between numerical analysis and analytical
% solution  by  \citet{deBoer1993}  in  terms  of a) solid displacements, b) fluid
% displacements,  c)  pore  pressure  and d) solid skeleton stress versus time for
% different depths.}}
% \end{center}
% \end{sffamily}
% \end{Large}
% \vspace{-.5cm}
% \end{figure}


% \begin{figure}[!htb]
% \begin{Large}
% \begin{sffamily}
% \begin{center}
% \includegraphics[width=1\textwidth]{FIGS/Hiremath_high.pdf}
% \caption{
% \label{solid_10_high} 
% Comparison   between   numerical   analysis   and   analytical   solution   
% \citet{Hiremath88}     for     the    case    of    strong    viscous    coupling
% ($k=0.148\times10^{-8}\text{~cm$^3$s/g}$)    in  terms of time histories of: a) solid velocity at
% $10$~cm  below  the  surface,  b) solid velocity at $30$~cm below the surface, c) fluid
% velocity  at  $10$~cm  below the surface and, d) fluid velocity at $30$cm below the
% surface.}
% \end{center}
% \end{sffamily}
% \end{Large}
% \end{figure}
% 
% \begin{figure}[!htb]
% \begin{Large}
% \begin{sffamily}
% \begin{center}
% \includegraphics[width=1\textwidth]{FIGS/Hiremath_low.pdf}
% \caption{\label{solid_30_high}{Comparison   between   numerical   analysis   and
% analytical  solution  \citet{Hiremath88}  for  the  case of weak viscous coupling
% ($k=0.148\times10^{-2}\text{~cm$^3$s/g}$):  a) solid
% velocity  at  $10$~cm  below  the  surface,  b) solid velocity at $30$~cm below the
% surface,  c)  fluid velocity at $10$~cm below the surface and d) fluid velocity at
% $30$cm below the surface, versus time.}}
% \end{center}
% \end{sffamily}
% \end{Large}
% \end{figure}
% 
% \clearpage
% 
% \begin{figure}[!htb]
% \begin{Large}
% \begin{sffamily}
% \begin{center}
% \includegraphics[width=1\textwidth]{FIGS/hiremath_NM06.pdf}
% \caption{\label{hiremath_NM06}{Analytical   solution  \citet{Hiremath88}  versus
% numerical  analysis,  with  two  different  combinations of spatial and temporal
% discretization  and numerical damping ($\gamma = 0.6$ and $\beta = 0.3025$), for
% the  case of weak viscous coupling ($k=0.148\times10^{-2}\text{~cm$^3$s/g}$): a)
% solid  velocity  at 10 cm below the surface and b) fluid velocity at 10 cm below
% the surface, with time.}}
% \end{center}
% \end{sffamily}
% \end{Large}
% \end{figure}

% 
% \begin{figure}[!htb]
% \begin{Large}
% \begin{sffamily}
% \begin{center}
% \includegraphics[width=1\textwidth]{FIGS/hiremath_fine.pdf}
% \caption{\label{hiremath_fine}{Analytical   solution  \citet{Hiremath88}  versus
% numerical  analysis  with  no  numerical  damping  ($\gamma  = 0.6$ and $\beta =
% 0.3025$),  in  terms  of fluid velocity time history at 10 cm below the surface,
% for the case of weak viscous coupling ($k=0.148\times10^{-2}\text{~cm$^3$s/g}$):
% a) $\Delta h = 0.001\text{~m}$ and $\Delta t = 2 \times 10^{-7}\text{~sec}$ and 
% b) $\Delta h = 0.005\text{~m}$ and $\Delta t = 5 \times 10^{-7}\text{~sec}$.}}
% \end{center}
% \end{sffamily}
% \end{Large}
% \end{figure}
% 
%        




\end{document}
\bye