\documentclass[review,authoryear,times]{elsarticle}

%\usepackage{lscape,graphicx}
%\usepackage{fancyheadings}
%\usepackage[numbers]{natbib}

%\setlength{\headrulewidth}{0.0mm}
%\setlength{\footrulewidth}{0.0mm}
%\lhead{}
%\rhead{}
%\chead{}
%\lfoot{}
%\cfoot{}
%\rfoot{\rm\thepage }
%%\addtolength{\headheight}{14pt}
%%\addtolength{\footheight}{14pt}
%
%\pagestyle{fancy}
%
%


\usepackage{url}
\newcommand{\ud}{{\rm d}}

%
%\def\baselinestretch{1.67}
%
%%12pt%12pt%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%12pt%12pt% this is for 12pt
%\setlength{\textwidth}{16.4cm}
%\setlength{\textheight}{22.8cm}
%\setlength{\hoffset}{-1.3cm}
%\setlength{\voffset}{-1.6cm}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%11pt%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%11pt% this is for 11pt
%%11pt\setlength{\textwidth}{15.2cm}
%%11pt\setlength{\textheight}{21.5cm}
%%11pt\setlength{\hoffset}{-0.5cm}
%%11pt\setlength{\voffset}{-0.5cm}
%

%\usepackage{times}
\usepackage{graphicx}
\usepackage{amsmath}
%\usepackage{color}
\usepackage{ulem}
%\usepackage{subeqn}
%


\begin{document}

\thispagestyle{empty}

\title{Numerical Modeling and Simulation of Pile in Liquefiable Soil}


\author[ZC]{
Zhao Cheng 
%(\raisebox{-1.0mm}
%{\includegraphics[width=3.55cm]{/home/jeremic/tex/works/IMENAuOriginaly/BorisJeremic.eps}}
}
%\ead{jeremic@ucdavis.edu}


\author[BJ]{
Boris Jeremi{\'c} 
%(\raisebox{-1.0mm}
%{\includegraphics[width=3.55cm]{/home/jeremic/tex/works/IMENAuOriginaly/BorisJeremic.eps}}
}
\ead{jeremic@ucdavis.edu}


%\author[rvt,focal]{K.?Bazargan\fnref{fn2}}
%\ead{kaveh@river-valley.com}

\cortext[cor1]{Corresponding Author: Boris Jeremi{\'c}}
%\cortext[cor2]{Principal corresponding author}
%\fntext[fn1]{This is the specimen author footnote.}
%\fntext[fn2]{Another author footnote, but a little more longer.}
%\fntext[fn3]{Yet another author footnote. Indeed, you can have
%  any number of author footnotes.}



%\address[focal]{
%Associate Professor, 
%Department  of  Civil  and  Environmental  Engineering,  University  of
%California, Davis, CA 95616.}
\address[ZC]{Staff Engineer, Earth Mechanics Inc. Oakland, CA 94621}

\address[BJ]{
Associate Professor, 
Department  of  Civil  and  Environmental  Engineering,  University  of
California, Davis, CA 95616.}


%
% \author{ 
% Zhao Cheng
% (\raisebox{-1.2mm}{\includegraphics[width=1.35cm]{/home/jeremic/tex/works/IMENAuOriginaly/Zhao.eps}})
% %(\raisebox{-1.2mm}{\includegraphics[width=1.1cm]{Zhao.eps}})
% \footnote{{Staff Engineer, Earth Mechanics Inc. Oakland, CA 94621}}
% %%%%%%%%
% \\
% and
% \\
% %%%%%%%%
% \setcounter{footnote}{3}
% Boris Jeremi{\'c} 
% (\raisebox{-1.0mm}{\includegraphics[width=3.55cm]{/home/jeremic/tex/works/IMENAuOriginaly/BorisJeremic.eps}})
% %(\raisebox{-1.0mm}{\includegraphics[width=3.2cm]{BorisJeremic.eps}})
% \footnote{Associate Professor, 
% Department  of  Civil  and  Environmental  Engineering,  University  of
% California, Davis, CA 95616. \texttt{jeremic@ucdavis.edu}, Coresponding Author}
% \\
% %Boris Jeremi{\'c} ({\cyr Boris Jeremi\cj{}})$^{1}$ \\
% %%%%%%%%
% }



\begin{abstract}

Presented  in  this  paper  is  numerical   methodology  to  model  and simulate
behavior of piles in liquefiable
soils.  Modeling  relies  on  use  of  validate
elasto-plastic  material  model for soil skeleton, verified fully coupled porous media
(soil  skeleton)  -- pore fluid (water) dynamic finite element formulation, and
detailed load staging of FEM models.  A  bounding  surface
elastic--plastic  sand  model that  accounts  for fabric change is used to model soil
skeleton, while a  fully coupled, dynamic, inelastic formulation  (u-p-U)  is used 
to model soil and water displacement and pore
water  pressures.
Much  attention is paid to accurate staged loading of the models, which start from
a  zero state of stress and strain for a soil without a pile, followed by application of
self  weight, then by excavation and pile installation with application of
pile  self-weighting.  Finally,  seismic  loading is applied followed by time to
dissipate  excess  pore pressures that have developed. A total of six cases were
modeled  and  simulated  varying slope inclination, presence of pile--column and
boundary condition for pile--column system.
%
Presented  are    interesting  and useful results that are used to  deepen our
understanding    of    behavior    of    soil--pile--column    systems    during
liquefaction (lateral deformations, pile pinning effect, ground settlement).  
%
Moreover,  detailed  description of modeling is used to emphasize
the  availability  and  use  of  high  fidelity  modeling  tools  for simulating
effects of liquefied soil on soil--structure systems.
\end{abstract}

%occurring,  excess pore water building-up and dissipation, soil skeleton lateral
%deformation,  ground  settlement,  pile response, and pile-soil interaction were
%interpreted at some length.

\begin{keyword} 
elasto-plastic, fully coupled finite  elements; liquefaction;  soil--pile dynamics; 
\end{keyword} 
 


\maketitle
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% \renewcommand{\baselinestretch}{1.4}
% \small\normalsize  % trick from Kopka and Daly p47.
% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
% \newpage
% \tableofcontents
% \newpage
% 
%%

%\newpage


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

\section{Introduction}  
\label{intro}






Liquefaction  is one of the most complex phenomena in
earthquake    engineering.  Liquefaction also represents  one of the biggest
contributors      to    damage  of  constructed  facilities  during  earthquakes
\citep{Kramer96}.
%
Prediction of  behaviors  of liquefiable soils  is difficult but achievable.
%
There  are  number  of  methods  that can be
utilized to predict such behaviors.
%
Methods currently used can have varying prediction accuracy and certainty.
%
%
Of  particular  interest  in this paper is the description of
verified and validated numerical simulation methodology based on
rational            mechanics      that   is used to model, simulate and predict
behaviors  of  a  single  pile  in liquefiable soil subjected to seismic
loading.
%
Both level and sloping ground pile systems are modeled and simulated.
%
Detailed    description   of  background  theory, formulation and implementation
were recently given by  \citet{Cheng2007} and \citet{Jeremic2007e}.



It  should be noted that presented development do show great promise in analyzing
a  myriad  of  liquefaction  related  problems  in  geotechnical  and structural
engineering. 
%
The  effectiveness  and  power  of  numerical  simulation  tools for analyzing
liquefaction  problems becomes even more important and prominent in the light of
potential disadvantages of models used in experimental simulations. 
%
These        disadvantages,   related  to  proper  scaling \citep{Wood2004} and 
problems  in  maintaining  appropriate similarities \citep{Harris1999} for first
order   important  phenomena,  can  render  scaled  models ineffective, when  used for physical
simulations (under one-step or multiple-step gravity loading).
%




In what follows, a brief literature review is provided.
%
The  literature  review comprises sections on observations of liquefaction
behavior    in   case  studies,   non--continuum modeling efforts,
review of 
redistribution    of    voids   and  pore  fluid  volume/pressures phenomena and
continuum
modeling efforts.
% %
% In  addition  to  that,  a  brief  comment on appropriateness of use of physical
% models  for  liquefaction  modeling  is  given  as well. 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\subsection{Background}  
%\label{background}
%
%
%
% %%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% \subsection{Literature Review}
% \label{LitRev}
% 
% 
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
% Behavior  of  piles  in  liquefiable  soil  has  been investigated for some time now.
% However,  numerical modeling and simulations of such behavior has been done only
% fairly recently.
% %
% It   is noted that inadequate performance of piles in liquefied soils
% has  been  the  main  cause  of  a number of bridge failures during earthquakes.
% %
% It  should be noted that, opposing to that, excellent  performance of pile
% foundation  systems  in  liquefied  soils  has  saved  quite  a few bridges (few
% examples noted below).
% 
% %

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\subsection{Observation of Behavior}
\paragraph{Observation of Behavior.}

Liquefaction  behavior  was observed during a number of earthquakes in the past.
%Listed below are selected references that describe behavior of  references that 
%
%
During  Alaskan  Earthquake (1964),  liquefaction was the main cause of
severe  damage  to  92  highway  bridges, moderate to light damage to another 49
highway  bridges,  and moderate to sever damage to 75 railroad bridges \citep{Youd89}.
%
During Niigata  Earthquake  (1964)  liquefaction induced
damage  to  foundation piles under Yachiya bridge \citep{Hamada92}.
%
During  that  same  earthquake,    girders    of  Showa  Bridge  toppled  as the
support   structure   and   piles   moved   excessively   due   to  liquefaction
\citep{JSCE66}.
%
During Kobe  Earthquake (1995), liquefaction was the primary cause of damage to
many  pile supported  or  caisson  supported  bridges  and  structures. 
%
For  example,  Shin--Shukugawa    bridge  was  subjected to excessive pile foundation
movement due to liquefaction \citep{Yokoyama97}.


Opposed    to  these  failures and collapses, there were a number of bridges with
pile  foundations  that  did  not suffer much or even minor damage even though there was
liquefaction around foundations.
%
For  example,  pile  foundations  of the Landing  Road  Bridge  in  New  Zealand
performed      quite    well   during  Edgecumbe   earthquake (1987) even with a
significant liquefaction recorded \citep{Berril97, Dobry2001}.
%
In  addition  to  that,  Second  Maya  Bridge piles (large steel pipes) were not
damaged  during  Kobe  earthquake despite significant liquefaction in surrounding
soils \citep{Yokoyama97}.







% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% \vspace*{-0.3cm}
% \paragraph{Piles in Liquefied Soils.}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\subsection{Non--Continuum Modeling Efforts}
\paragraph{Non--Continuum Modeling Efforts.}

Modeling    and   simulation  of  piles in liquefied grounds has been focus of a
number of recent studies. 
%
The  simple approach, based on scaling of p-y springs has been suggested early
by   \citet{JRA80},  \citet{AIJ88},  \citet{Liu95},  \citet{Miura89}    and
\citet{Miura91}.
%
However, large inconsistencies with material parameter selection are present when
p-y spring approach is used for piles in liquefied soils. 
%
Since  p-y  methodology for liquefied soils is not based on rational mechanics,
appropriate choice of material parameters is primarily based on
empirical  observations  of behaviors of piles in liquefied soils in experimental
studies. 
%
A  number  of experimental  studies have carefully examined pile behaviors in
liquefiable soils. We mention  \citet{Tokida92}, \citet{Liu95}, \citet{Abdoun97},
\citet{Horikoshi98} and \citet{Boulanger06}.
%
Studies using physical model can be used to obtain very high quality data on behavior of piles
in  liquefied soils, provided that similarity of important physical phenomena is
maintained \citep{Wood2004, Harris1999}. 
%
% %Comments  on physical modeling of models is provided in some more
% %details below in section \ref{Comments}. 
% 
% 
% 
% 
% In recent years, numerical modeling of piles in liquefied soils has 
% relied on beam models on soil springs.
Some   of  the  recent  papers  that  discussed  use  of  these  models  and  gave
recommendations about parameter choices are listed for reference:
\citet{Tokimatsu98}, 
\citet{Martin02}, 
\citet{Dobry03}, 
\citet{Liyanapathirana05a},
\citet{Rollins05}, 
\citet{Cubrinovski06},
\citet{Brandenberg07a}.  
%\citet{Tokimatsu03}, 


% %Some  of  the differences in
% %recommended  parameter  relationships  and  loading  conditions are important to
% %practice and warrant further study to resolve.
% 
% The   currently proposed   simulation (design)   guidelines   are   mainly  for 
% two major approaches \citep{Boulanger07}:  (1) piles  in laterally spreading ground where the
% out-of-plane    thickness    of  the  spreading  ground  is  sufficiently large
% that free-field  soil  displacements are relatively unaffected by the presence of the
% piles and (2) piles at approach embankments where the restraining or
% \emph{pinning}  effects  of  the  piles and bridge superstructure can reduce the
% lateral  spreading  demands imposed by the finite width embankment. Guidance for
% both situations is limited to equivalent static analysis using beam on nonlinear
% Winkler  foundation  (BNWF)  analysis  methods,  with recognition that nonlinear
% dynamic analysis may be warranted for important bridges.
% 
% The  static  BNWF  analysis  methods  has  been proposed in two major forms: (a)
% BNWF--SD  with  imposed  soil  displacements and (b) BNWF--LP with imposed limit
% pressures,    The        first  approach is preferred   by some researchers, most
% notably by
% \citet{Boulanger07}.  In  reality,  however,  for both of these cases there is a
% coupling  between  the  pile  structure and the surrounding soil which in result
% determines  the  so  called  imposed  soil  displacements  and  limit pressures.
% %Moreover,   there   are   significant   uncertainties   involved  in  predicting
% %liquefaction--induced ground displacements for lateral spreads.
% 

% 
% While  the  problem  of  piles in liquefiable soils has only recently started to
% receive  increased  attention,  there  is already sizable body of work, as noted
% above.  
% %
% %What  is  missing, and is one of the main parts of this proposal, is
% %the  development  of  methodology that will allow for simulations of problems of
% %piles  in  liquefied  soils,  using (already existing) state  of  the  art
% %computational geomechanics formulations and implementation. 
% %
% %
% %%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\subsection{Redistribution of Voids and Pore Fluid Volume/Pressures}
\paragraph{Redistribution of Voids and Pore Fluid Volume/Pressures.}



Mechanics  of pile behavior in
liquefiable grounds is based on the concept of redistribution    of  voids  and 
pore fluid volume/pressures (RVPFVP).
%
It should be  emphasized that geomechanics phenomena of redistribution
of voids -- pore fluid
volume/pressure is used here in purely mechanistic way. 
%
That is,  RVPFVP is a phenomena
that occurs in saturated soils and that phenomena is responsible for (is manifested in)
liquefaction  related  soil behaviors with or without piles.
%
%
This is noted as in some recent publications,  RVPFVP terminology
is explicitly used  for problems
of liquefaction induced failures of sloping grounds without piles.
%
Our understanding of the  RVPFVP  phenomena is that RVFVP is responsible for
many    more    facets   of behavior of liquefied soils, 
rather than only failure of liquefied slopes.



The  early  investigation  of  the  RVPFVP  phenomena  was related to the
behavior of infinite slopes.
%
For example,  loss of shear strength in infinite slopes is one of the
early understood manifestations of RVPFVP  \citep{Whitman85,NRC85,Malvick06}.
%
Laboratory  investigation  of sand was also used to observe the RVPFVP phenomena
\citep{Casagrande78,Gilbert84}



% In this particular case  of liquefiable slopes, parts of liquefied soil will locally loosen
% while  other   parts will densifies  \citep{Whitman85}, \citep{NRC85}. The loosening and
% subsequent strength  loss  (from RVPFVP) is a consequence of  water
% seepage  due  to  development of excess  pore  pressure  and  restriction of flow due to
% local  permeability  conditions.  The slope can undergo significant deformations
% if the loss of strength associated with loosening is large enough to reduce the
% undrained strength below levels required for stability of the
% slope \citep{Malvick06}.









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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\subsection{Continuum Modeling Efforts}
\paragraph{Continuum Modeling Efforts.}

Continuum  based formulations for  modeling liquefaction problems have been
present  for over two decades. 
%
In   a   landmark   paper,   \citet{Zienkiewicz84}   presented  three  possible
coupled formulations for modeling of soil skeleton -- pore fluid problems.
%
The  most  general  and complete one is the so called u-p-U formulations while
the  other  two,  the  u-p  and the u-U have a number of restrictions on the
domain of application. Here, the unknowns are  the soil skeleton displacements
$u$;    the  pore  fluid  (water)    pressure  $p$;   and the pore fluid (water)
displacements $U$. The u-p formulation captures the
movements of the soil skeleton and the change of the pore pressure, and is the most
simplistic  one of the three mentioned above. This formulation neglects the
differential  accelerations of the pore fluid (it does account for  acceleration
of  pore  fluid  together  with  soil  skeleton,  but not the separate one if it
exists), and in one version neglects the compressibility
of  the  fluid  (assuming  complete incompressibility of the pore fluid). In the
case of incompressible pore fluid, the formulation requires
special  treatment  of  the  approximation  function (shape function)  for pore fluid to
prevent the volumetric locking  \citep{Ziekienwicz2000a}. 
The  majority  of  the  currently  available  implementations  are based on this
formulation. For example \cite{Elgamal2002} and \cite{Elgamal2003}
developed  an  implementation  of  the  u-p formulation with the multi-surface
plasticity model by \cite{Prevost85}, while \cite{Chan88}  and \cite{Zienkiewicz98}
used  generalized theory of plasticity \cite{Pastor90}.
%

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 variables, 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 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 additional 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. With known (given)
volumetric  compressibility  of  the  pore fluid, pore fluid pressure can be calculated.
%
Despite  it's  power,  the  u-p-U formulation has rarely been implemented into
finite  element  code,  and  has  never (at least to our knowledge) been used to
analyze pile -- liquefied soil interaction.
%
This  can  be  attributed  in  part  to  a  sophistication of  implementation that is
required, and to a sizable increase in computational cost for u-p-U elements.



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Numerical Formulation, Verification and Validation}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%




\subsection{Finite Element Formulation }




The  discretized,  finite element system of equations for u-p-U formulation,
which is   based   on   earlier   work   by \citet{Zienkiewicz84}, 
can be written
%as
%%%
%\begin{eqnarray}
%\left[ \begin{array}{ccc}
%\bf{M_s} & 0 & 0 \\
% 0  & 0 & 0 \\
% 0  & 0 & \bf{M_f}
%\end{array} \right]
%\left[ \begin{array}{c}
%\bf{\ddot{\bar{u}}} \\
%\bf{\ddot{\bar{p}}} \\
%\bf{\ddot{\bar{U}}}
%\end{array} \right]
%+
%\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]
%\left[ \begin{array}{c}
%\bf{\dot{\bar{u}}} \\
%\bf{\dot{\bar{p}}} \\
%\bf{\dot{\bar{U}}}
%\end{array} \right] \nonumber\\
%+
%\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]
%\left[ \begin{array}{c}
%\bf{\bar{u}} \\
%\bf{\bar{p}} \\
%\bf{\bar{U}}
%\end{array} \right]
%+
%\left[ \begin{array}{c}
%\bf{f_{in}^u(u)} \\
%0  \\
%0
%\end{array} \right]
%-
%\left[ \begin{array}{c}
%\bf{\bar{f}^u } \\
%0  \\
%\bf{\bar{f}^U}
%\end{array} \right]
%=
%0
%\label{67}
%\end{eqnarray}                                                                                                                             %%
%or
in tensor index form
%%   %%
%%   \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 \Omega \\
%%   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}
%%   %%
%%
%%
%%
\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] 
+
\left[ \begin{array}{c}
\int_{\Omega} N_{K,j}^u \sigma''_{ij} \ud \Omega \\
0  \\
0
\end{array} \right]
-
\left[ \begin{array}{c}
\bar{f}^u_{Ki} \\
0  \\
\bar{f}^U_{Ki}
\end{array} \right]
= 0
 \nonumber \\ 
\label{68}
\end{eqnarray}
%%
%%where
%%%%
%%\begin{equation}
%% \bf{ f_{in}^u(u)} = \int_{\Omega} N_{K,j}^u \sigma''_{ij} \ud \Omega
%%\end{equation}                                                                                                                           %%
where the components are
%%%     %%%
%%%     \begin{eqnarray}
%%%     % \bf{M_s} &=&
%%%      (M_s)_{KijL} &=& \int_{\Omega} N_K^u (1-n) \rho_s \delta_{ij} N_L^u \ud \Omega \nonumber\\
%%%     % \bf{M_f} &=&
%%%      (M_f)_{KijL} &=& \int_{\Omega} N_K^U n \rho_f \delta_{ij} N_L^U \ud \Omega \nonumber\\
%%%     % \bf{C_1} &=&
%%%      (C_1)_{KijL} &=& \int_{\Omega} N_K^u n^2 k_{ij}^{-1} N_L^u \ud \Omega \nonumber\\
%%%     % \bf{C_2} &=&
%%%      (C_2)_{KijL} &=& \int_{\Omega} N_K^u n^2 k_{ij}^{-1} N_L^U \ud \Omega \nonumber\\
%%%     % \bf{C_3} &=&
%%%      (C_3)_{KijL} &=& \int_{\Omega} N_K^U n^2 k_{ij}^{-1} N_L^U \ud \Omega \nonumber\\
%%%     % \bf{G_1} &=&
%%%      (G_1)_{KiM} &=& \int_{\Omega} N_{K,i}^u (\alpha-n) N_M^p \ud \Omega \nonumber\\
%%%     % \bf{G_2} &=&
%%%      (G_2)_{KiM} &=& \int_{\Omega} n N_{K,i}^U N_M^p \ud \Omega \nonumber\\
%%%     % \bf{P}   &=&
%%%      P_{NM} &=& \int_{\Omega} N_N^p \frac{1}{Q} N_M^p \ud \Omega
%%%     \label{69-1}
%%%     \end{eqnarray}
\begin{eqnarray}
% \bf{M_s} &=&
 (M_s)_{KijL} = \int_{\Omega} N_K^u (1-n) \rho_s \delta_{ij} N_L^u \ud \Omega
 \;\;&\mbox{;}&\;\;
% \bf{M_f} &=&
 (M_f)_{KijL} = \int_{\Omega} N_K^U n \rho_f \delta_{ij} N_L^U \ud \Omega \nonumber\\
% \bf{C_1} &=&
 (C_1)_{KijL} = \int_{\Omega} N_K^u n^2 k_{ij}^{-1} N_L^u \ud \Omega
 \;\;&\mbox{;}&\;\;
% \bf{C_2} &=&
 (C_2)_{KijL} = \int_{\Omega} N_K^u n^2 k_{ij}^{-1} N_L^U \ud \Omega \nonumber\\
% \bf{C_3} &=&
 (C_3)_{KijL} = \int_{\Omega} N_K^U n^2 k_{ij}^{-1} N_L^U \ud \Omega
 \;\;&\mbox{;}&\;\;
% \bf{G_1} &=&
 (G_1)_{KiM} = \int_{\Omega} N_{K,i}^u (\alpha-n) N_M^p \ud \Omega \nonumber\\
% \bf{G_2} &=&
 (G_2)_{KiM} = \int_{\Omega} n N_{K,i}^U N_M^p \ud \Omega 
% \bf{P}   &=&
 \;\;&\mbox{;}&\;\;
 P_{NM} = \int_{\Omega} N_N^p \frac{1}{Q} N_M^p \ud \Omega \nonumber \\
% \bf{f^u} &=&
 (\bar{f}_s)_{Ki} = (f_1^u)_{Ki} - (f_4^u)_{Ki} + (f_5^u)_{Ki}
% \bf{f^U} &=&
 \;\;&\mbox{;}&\;\;
 (\bar{f}_f)_{Ki} = - (f_1)_{Ki} + (f_2)_{Ki}  \nonumber\\
%
 (f_1^u)_{Ki} =
 \int_{\Gamma_t} N_K^u \sigma''_{ij} r_j \ud \Gamma 
 \;\;&\mbox{;}&\;\;
 (f_4^u)_{Ki} =
 \int_{\Gamma_p} N_K^u (\alpha-n) p r_i \ud \Gamma \nonumber\\
%
 (f_5^u)_{Ki} =
 \int_{\Omega} N_K^u (1-n) \rho_s b_i \ud \Omega 
 \;\;&\mbox{;}&\;\;
 (f_1^U)_{Ki} =
 \int_{\Gamma_p} N_K^U n p r_i \ud \Gamma \nonumber\\
%
 (f_2^U)_{Ki} =
 \int_{\Omega} N_K^U n \rho_f b_i \ud \Omega
\label{69-1}
\end{eqnarray}
%%
Here  $N^u$,  $N^p$,  $N^U$ are shape functions of solid skeleton displacement, pore pressure
and fluid displacement, respectively;
$\rho$, $\rho_s$, $\rho_f$ are density of the total, solid and fluid phases, respectively;
$n$  is  porosity,
the symbol $r_i$ is the direction of the normal vector on the boundary;
$\bar{u}_{Lj}$ are the nodal displacements of the solid part;
$\bar{p}_M$ are the nodal  pore pressures and
$\bar{U}_{Lj}$ are the nodal displacements of the fluid part.
%
Moreover, $\Omega$ represents the domain of interest; $\Gamma_t$ is the traction boundary,
and $\Gamma_p$ is the pressure boundary.
%
%
Set    of  dynamic  Equations (\ref{68}) represents the most general formulation
and discretization for a material nonlinear (inelastic)  porous medium (soil skeleton)  
that  is  fully  saturated  with  linear elastic, compressible pore fluid
(water).
%
%There  are number of feature of this formulation and discretization that deserve
%further comments.
%
%


Water    accelerations    are    explicitly    taken   into  account,  both  for
conforming  and   differential movements with respect to the soil skeleton. This
proves  to  be  very  important  for  models  where  porous  soil is adjacent to
structural foundations (piles for example). In these models, the dynamics of two
model components, saturated soil and piles, are quite different and there exists
a possibility of significant relative movement (displacements, velocities and
accelerations) between soil skeleton and the concrete piles or footings. If this
relative movement exists, pore water is pumped in and out of soil skeleton, thus
creating   a  significant  differential  accelerations  relative  to  the  soils
skeleton. 

Both   velocity  proportional  and  displacement  proportional  damping  follow
directly from formulation and discretization. 
%
The  velocity  proportional damping is taken into account through damping matrix
in Equation (\ref{68}). Of particular interest are sub-matrices (tensors) 
$(C_1)_{KijL}$,   $(C_2)_{KijL}$   and  $(C_3)_{KijL}$  (Equations~(\ref{69-1})).
Physically,  those  sub-matrices  represent  coupling  of  pore  water  and solid
skeleton,  which  is  velocity  proportional.  This  coupling  is  a function of
permeability  $k$  and  porosity  $n$. For example, from Equation~(\ref{69-1}), it
follows   that  for a  soil with larger  porosity  $n$, 
which has more  pores  in a soil skeleton and therefore
more  pathways  for  pore  water to travel through the soil, the damping will be
higher.  Similarly,  for  a soil with smaller permeability, where pore water has
more  difficulty  in  traveling  through  pores,  where there is more friction of
flowing  pore  fluid  with  the  soil skeleton, the velocity proportional
damping will be increased.
%
Displacement  proportional damping is controlled by inelastic
material   behavior,   which  in  turn  is  controlled  in  the  finite  element
discretization  by  the  nonlinear,  elastic--plastic stiffness matrix, given as
$\int_{\Omega}      N_{K,j}^u      \sigma''_{ij}      \ud   \Omega$      (see   
Equation~(\ref{69-1})) here in order to correctly
account for general, nonlinear dynamic time integration \citep{local-87}.
 





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


\subsubsection{Time Integration}
\label{Time Integration}


Numerical    integration   of  Equations  of motion (Equation \ref{68}) is done
using     Newmark        \citep{Newmark59}        or    Hilber-Hughes-Taylor
\citep{HHT,Hughes1978a,Hughes1978b}  algorithms. Of particular importance is the
proper  algorithmic  treatment  for  nonlinear analysis which
introduces  changes  in  a  way  residual forces are calculated \citep{local-87,
Jeremic2007e}.

The finite element discretization in Equation~(\ref{68}) defines a damping matrix,
which takes into account physics of velocity dependent interaction of pore water
and   soil skeleton. This damping matrix is more appropriately used here than damping
matrix introduced through Rayleigh damping \citep{Chopra01}.
It  should  be  note  that a small amount of numerical damping is used with both
Newmark  and  Hilber-Hughes-Taylor methods, in order to damp out response in higher frequencies 
that    is   introduced   by   the   spatial   finite   element   discretization
\citep{local-86}.





% 
%  
%  
%     Equation \ref{68} can be written in a simplified residual form
%     %%%
%     \begin{eqnarray}
%       R = M \ddot{x} + C \dot{x} + K' x + F(x) - f = 0
%     \label{71}
%     \end{eqnarray}
%     %%%
%     where $x = \{\bar{u}_{Lj}, \; \bar{p}_M, \; \bar{u}_{Lj}\}^T$ represents the generalized unknown variables.
%     Equation \ref{71} represents a general non-linear form.
%  
%     Numerical damping introduced in the Newmark time integration method
%     will degrade the order of accuracy.
%     The Hilber-Hughes-Tailor (HHT) time integration $\alpha$-method
%     (Hilber et al \cite{HHT}, Hughes and Liu \cite{Hughes1978a} and \cite{Hughes1978b})
%     using an alternative residual form
%     by introducing an addition parameter $\alpha$ to improve the performance.
%     The current residual at time $n+1$ (for simplicity, the left superscripts $n+1$ are omitted)
%     can be estimated from the previous variables at time $n$ by:
%     %%%%
%     \begin{eqnarray}
%       R = M \; \ddot{x}
%                  + (1 + \alpha) [C \dot{x} + K' x + F(x)]
%                  - \alpha [C\;{}^{n}\dot{x} + K'\;{}^{n}x + F({}^{n}x)]
%                  - \; f
%     \label{upU_HHT_R}
%     \end{eqnarray}
%     %%%%
%     but retaining the Newmark finite-difference formulas
%     and parameters $\beta$ and $\gamma$:
%     %%
%     \begin{eqnarray}
%       x &=& {}^n x + \Delta t \; {}^n \dot{x} + \Delta^2 t \; [(\frac{1}{2} - \beta) \; {}^n \ddot{x} + \beta \; \ddot{x}]
%        \label{upU_Newmark_EQ1} \\
%       \dot{x} &=& {}^n \dot{x} + \Delta t \; [(1 - \gamma) \; {}^n \ddot{x} + \gamma \; \ddot{x}]
%       \label{upU_Newmark_EQ2}
%     \end{eqnarray}
%     %%
%     If $\alpha = 0$,
%     Equation \ref{upU_HHT_R} is reduced to
%     the  special case of the HHT time integration method
%     which is exactly the Newmark time integration method.
%     Decreasing $\alpha$ value increase numerical dissipation (Hughes \cite{local-86}).
%  
%     The iteration procedure is
%     \begin{eqnarray}
%       \left[ M + (1+\alpha) \gamma \Delta t C + (1+\alpha) \beta \Delta t^2 K \right] \Delta \ddot{x} = - R
%       , \;\;\;
%       \ddot{x} = \; {}^n \ddot{x} + \Delta \ddot{x}
%     \label{upU_HHT_Newton}
%     \end{eqnarray}
%     %%%
%     where
%     \begin{eqnarray}
%       K = K' + \frac{\partial F}{\partial x}
%     \end{eqnarray}
%  
%     If the parameters $\alpha$, $\beta$ and $\gamma$ satisfy
%     %%%%
%     \begin{equation}
%       -\frac{1}{3} \le \alpha \le 0, \;\;\; \gamma = \frac{1}{2}(1-2\alpha), \;\;\; \beta = \frac{1}{4}(1-\alpha)^2
%     \label{upU_HHT_ab}
%     \end{equation}
%     %%%
%     it is unconditionally stable and second-order accurate (Hughes \cite{local-86}).
%  
%  %   There are different denotation meaning of the parameter $\alpha$,
%  %   e.g. the parameter $\alpha$ in the HHT integrator codes of OpenSees
%  %   equals to the conventional $\alpha$ plus one.
%  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{Material Model}
\label{Material_Model}


Material model plays one of the key roles in numerical simulation of the dynamic
behaviors  of liquefiable soil. Correct modeling of volumetric response by a good
model  on  the  constitutive  level allows for accurate modeling of the boundary
value problems where liquefaction is involved. 
%
In this work, a critical state soil mechanics based model developed by
\citet{Manzari97} and \citet{Dafalias04} is used. 
%
%This model considered the effects of the state parameter on the behaviors
%of the dense or loose sands.
Among    many    excellent    features   of this model we note the capability to
utilize  a  single      set    of material   parameters for a wide range of void
ratios and stress states for the same soil.  This  feature  allows  
the same  material parameters to be used from the very
beginning    of     loading  (from  zero  stress/strain  state),  through  self 
weight,  pile  construction  and  finally  dynamic shaking. In addition to that,
model  validation  for  Toyoura  sand,  used  in  this  study, shows excellent
agreement with test data \citep{Jeremic2007e}.
%
It is emphasized again that  Dafalias-Manzari set of models used here, together 
with  powerful  constitutive  integrations,
u-p-U   formulation   and  proper  dynamic  equation  of  motion time integration
\citep{Jeremic2007e},      is    capable  to  follow  soil  material  from  zero
stress/strain  state      (before    applying  self  weight)  all  the  way  to 
pre--liquefaction, liquefaction and post liquefaction behaviors. This transition
%%
\begin{itemize}
%
\item       from      pre--liquefaction,      where      mean  effective  stress
$p^{\prime}(=-\sigma^{\prime}_{kk}/3)$ is being continuously reduced (as pore water
pressure is increased) to become numerically almost zero,
%
\item  through  liquefaction, where mean effective stress $p^{\prime}$  is  numerically 
almost zero,
%
\item  to  post  liquefaction,  where  mean  effective  stress  $p^{\prime}$  is
increasing (as pore water pressure is decreasing due to diffusion),
%
\end{itemize}
%%
is  successfully  handled  by  the  material  model  with  the  help  of proper
constitutive  and  dynamic level finite element integrations. Each of the above
phases  results  in  proper change (evolution) of the void
ratio and the soil skeleton fabric, representing material internal variables.





%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Verification and Validation }
\label{VandV}%


Confidence  in  numerical  prediction  is  firmly based on a valid
process of verification and validation \citep{Oberkampf2002}. 
%
Verification provides
evidences  that  the  model  is  solved  correctly. Verification is also meant to
identify    and remove errors in computer coding and verify numerical algorithms
and is desirable in quantifying numerical errors in computed solution.
%
Validation provides evidences
that    the    correct  model  is  solved. It is custom to identify verification
with mathematics of the problem and validation  with the physics of the problem.

Computational    simulation  tools      used  in  this  study underwent detailed
verification and validation program.
% 
Validation  consisted of applying a comprehensive   set of software verification
tools,  available  in  public  domain  and/or  developed  at  the  Computational
Geomechanics Lab at UC Davis to developed libraries and programs.
%%for removing errors in implementation. 
%
In  addition to that, a number of closed form solutions were used to verify that
our  models were solved correctly. A set of problems with solutions
developed      by      \citet{Coussy95}      were    numerically    modeled  and
simulated during validation process \citep{Jeremic2007e}.
%
Mentioned are verification problems of 1D consolidation, 
line injection of fluid in a reservoir, and cavity expansion in saturated medium
(2D and 3D).

In  addition  to  those  problems  which  dealt  with  relatively slow phenomena
without  significant  influence of inertial forces, a truly dynamic problem of
shock  wave  propagation  in  porous  medium  was  also  used  for verification 
\citep{Gajo1995,Gajo1995b}.
%
Verification  examples  using  shock  wave propagation are most
sever  tests of the u-p-U formulation and the numerical integration
algorithms for the dynamic equations of motions.
%
It    should   be  noted that  these verification test provided excellent, close
matching  of numerical and closed for solutions, within the known limitations
of numerical accuracy of finite element discretization.
%
This  limitation  in  accuracy  is  a      simple and expected  consequence  of 
approximate  nature of the finite element method \citep{FEM-IV-1, FEM-IV-2}, and
the    finite    precision    arithmetic    used    in   computer   calculations
\citep{Schnabel83}.



Validation was done by simulating constitutive behaviors of sand material used in
predictions   (Toyoura sand). Comparison of validation predictions with physical
test   data   \citep{Jeremic2007e}   shows   very  good  predictive
capabilities  of  the  material  model  as  well  as of the underlying numerical
integration algorithms on the constitutive level. 




\section{Staged Simulation Model Development}
%\section{Simulation Model Development}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

Model  development for a pile in the liquefiable soil follows physics (mechanics) of
the problem as close as possible.
%
Numerical  simulation  of  such problems in geomechanics is usually based on stages of
loading and increments within those stages.
% 
%Hence, modeling follows the load staging in detail.
%
% Load  stages  are  centered around the assumed past history of loading on a soil
% pile system. 
% 
% 
% 
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% \subsection{Staged Model Construction}





All  load  stages  are  applied  to  a  series   of finite element models, all of
which share features of  an initial soil model.
% 
This  initial  soil  model  consists of a soil block with dimension  of
$12 \times 12 \times 15$~m (length $\times$ width $\times$ depth).
%
Due  to  the  symmetry  of  the model, only half of the block is modeled. 
%
Symmetry    assumptions    is based on assumption that  all  the  loads, dynamic
shaking and other influences are symmetric with respect to the plane of symmetry.
%
This specialization
to  symmetric  model  reduces  model  generality (for example this use of symmetry will
preclude  analysis  of  dynamic  shaking  perpendicular  to sloping ground dip).
%
However, as   our  goal is to present a methodology of analyzing behavior of piles in
liquefying  ground,  this  potential  drawback  is  not deemed significant in this
study.
%
Finite element mesh for the model is presented in Figure~(\ref{FEmesh3D}).
The  initial mesh  consists  of  160 eight node u-p-U elements. 
%


%
\begin{figure}[!htbp]
\begin{center}
\includegraphics[width=0.43\textwidth]{NewFiga/FEmesh3D.eps}
\hfill
\includegraphics[width=0.52\textwidth]{NewFiga/FEmeshPileBeam.eps}
\caption{
  \label{FEmesh3D}  
  {Left: Three dimensional finite element mesh featuring initial soil setup, where all
  the  soil  elements  are  present.  The  gray region of elements is excavated
  (numerically)  and  replaced  by a pile during later stages of loading; Right:
  Side  view  of the pile-soil model with some element and node annotation, used
  to visualize results.}
}  
\end{center}
\end{figure}
%


Each  node  of  the  mesh  has  7  degrees of freedom, three for soil
skeleton displacements  ($u_i$),  one for pore water pressure ($p$), and three
for pore water displacement ($U_i$).
%
While it can be argued that  the  mesh  is  somewhat  coarse, it is well refined
around the pile, yet to be installed, in place of gray region in the middle.
%

%%   
\begin{table}[tbh]
\begin{center}
\caption{ Material parameters used for Dafalias-Manzari elastic--plastic model. }
% {\normalsize
\begin{tabular}
[c]{|l|c|c||l|c|c|}\hline
\multicolumn{2}{|l|}{Material Parameter} & Value &
\multicolumn{2}{|l|}{Material Parameter} & Value\\\hline
Elasticity        & $G_{0}$         & 125~kPa    & Plastic modulus  & $h_{0}$    & 7.05\\     
                  & $v$             & 0.05       &                  & $c_{h}$    & 0.968\\    
Critical sate     & $M$             & 1.25       &                  & $n_{b}$    & 1.1\\      
                  & $c$             & 0.8        & Dilatancy        & $A_{0}$    & 0.704\\    
                  & $\lambda_{c}$   & 0.019      &                  & $n_{d}$    & 3.5\\      
                  & $\xi$           & 0.7        & Fabric-dilatancy & $z_{max}$  & 4.0\\
                  & $e_{r}$         & 0.934      &                  & $c_{z}$    & 600.0\\
Yield surface     & $m$             & 0.02       &                  &            & \\\hline
\end{tabular}
% }
\end{center}
\label{DMTest_parameters}%
\end{table}
%%   
%%
   
%%   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{table}[tbh]
\begin{center}
\caption{ \label{OneC_Material_Parameter}
Additional parameters  used  in  FEM simulations.}
\begin{tabular}[c]{|c|c|c|} \hline
\multicolumn{2}{|c|}{Parameter} & Value \\  \hline 
Solid density                    & $\rho_s$      & $2800~kg/m^3$ \\
Fluid density                    & $\rho_f$      & $1000~kg/m^3$ \\
Solid particle bulk modulus      & $K_s$         & $1.0 \times 10^{12}~kN/m^2$ \\
Pore fluid bulk modulus               & $K_f$         & $2.2 \times 10^{6}~kN/m^2$ \\
%Biot coefficient                 & $\alpha$      & $1.0$  \\
permeability                     & $k$           & $1.0\times 10^{-4}$~$m/s$           \\
Gravity                          & $g$           &  10~m/s$^2$ \\ \hline
%
\end{tabular}
\end{center}
\end{table}
%  

A  single  set  of parameters  is used with the Dafalias-Manzari material model.
Soil  is  modeled  as  Toyoura sand and material parameters (summarized in Table
\ref{DMTest_parameters})     are     calibrated     using     tests    by
\citet{Verdugo1996}, while initial void ration was set to $e_0 = 0.80$.
%
It  is  very  important to emphasize that the state of stress
and  internal  variables from initial state (zero for stress and given value for
void ratio and fabric) will evolve through all stages of loading 
by proper modeling and algorithms, 
by using single set of material parameters.
%
%
%
Table  \ref{OneC_Material_Parameter}  presents additional parameters,  other  than material
parameters  presented  in  Table  \ref{DMTest_parameters},  used  for  numerical
simulations.



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{First Loading Stage: Self Weight}


The initial stage of loading is represented by the application of self weight on
soil, including both the soil skeleton and the pore water.
%
Initial state in soil before application of self
weight is of a zero stress and
strain  while  void ratio and fabric are given initial
values. 
%
The state of stress/strain, void ratio and fabric will evolve upon application of self weight.
%
At  the  end  of  self  weight loading stage, soil is under appropriate state of
stress  ($K_0$ stress),  the  void  ratio corresponds to the void ratio after self
weight (redistributed such that soil is denser at lower layers),
while soil fabric has evolved with respect to stress induced anisotropy.
%
All  of  these  changes  are  modeled  using Dafalias--Manzari
material  model  and  using  constitutive  and  finite element level integration
algorithms  developed within UC Davis Computational Geomechanics group in recent
years.
%

Boundary  conditions (BC) for self weight stage of loading are set in the following
way:
%

\begin{itemize}



\item  Soil  skeleton  displacements ($u_i$), are fixed  in all three directions
at    the    bottom   of the model.
At the side planes, nodes move only vertically to mimic self-weight effect. 
All other nodes are free to move in any direction. 


%
%free to move only vertically at four
%sides  of  the  model, and are free to move in all three directions at all other
%nodes.


\item      Pore  water  pressures ($p$), are free to develop at the bottom plane and at
all    levels    of  the  models except at the top level  at soil 
surface where they are fixed (set to zero, replicating drained condition),


\item   Pore water  displacements ($U_i$), are fixed  in all three directions at
the bottom, are free to move only vertically at four sides of the
model   and are free to move in any direction at all other nodes.



\end{itemize}
%
%
These  boundary condition are consistent with initial self-weighting deformation 
condition for soil and pore water at the site.


For   the  case of sloping ground, an additional load sub-stage is applied after self
weight loading, in order to mimic self weight of inclined (sloping) ground.
%
This  is effectively achieved by applying  a  resultant  of total self weight of
the  soil  skeleton  times the sine of the inclination angle at uphill side of
the  model. This load is applied only to the solid skeleton DOFs, and not on the
water  DOFs.  Physically  it  would  be  correct  to consider the sloping ground
effects  on  the  pore water as well. This will create a constant flow field of
the  water downstream, which, while physically accurate, is small enough that it
does not have any real effect on modeling and simulations performed here.
%



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Second Loading Stage: Pile--Column Installation}
\label{pile_installation}




After  the first  loading stage,  comprising  self  weight applications (for level or sloping
ground,    as    discussed    above), second loading stage includes installation
(construction) of the pile--column.
Modeling changes performed during loading stage included:

\begin{itemize}

\item  Excavation of soil occupying space where the pile will be installed. This
was done by removing elements, nodes and loads on elements shown
in gray in Figure~(\ref{FEmesh3D}).

\item    These    elements   were  replaced  by very soft set of 
elements with small
stiffness,  low permeability. This was done in order to prevent water from
rushing  into  the  newly opened hole in the ground after original soil elements
(used in the first loading stage) are removed. 


\item    Installation of a pile in the ground and a superstructure (column) above the 
ground. Nonlinear  bean--column elements were used for both pile and column
together    with  addition of appropriate nodal masses at each beam-column node,
and  with  the  addition of a larger mass at the top representing lumped
mass  of  a bridge superstructure. Pile beam-column elements were connected with
soil skeleton part of soil elements using a specially devised technique.


\end{itemize}


As  mentioned  earlier, the volume that would be physically occupied by the pile
in the pile hole,   
is ``excavated" during  this loading stage. 
%
Beam--column elements, representing piles, are then placed in the middle of this
opening.  Pile (beam--column)  elements   are then connected  to  the
surrounding  soil elements by means of stiff elastic beam--column
elements.  These ``connection" beam--column elements extend from each pile
node   to  surrounding  nodes  of  soil elements.  The
connectivity      of  nodes to soil
skeleton nodes is  done  only for three beam--column translational DOFs,
while  the  three rotational DOFs from the
beam--column   element   are  left  unconnected.  These  three  DOFs  from  the
beam--column  side  are  connected  to  first  three  DOFs  of  the u-p-U soil
elements, representing displacements of the soil skeleton ($u_i$).
%
Water  displacements ($U_i$) and pore water pressures ($p$) are not connected in
any  way.  Rather,  these  two  sets of DOFs representing pore water behave in a
physical  manner  (cannot  enter  newly  created  hole  around pile beam--column
elements)  because  of the addition of a soft, but very impermeable set of u-p-U
elements,  replacing  excavated  soil elements. By using this method, both solid
phase (pile with soil skeleton) and the water phase (pore water within
the soil) are appropriately modeled.
%
%
Figure~(\ref{coupling})      shows in some detail schematics of coupling between
the pile and soil skeleton part finite elements.

%                          
%
\begin{figure}[!htbp]
\begin{center}
%{\includegraphics[width=8cm]{/home/jeremic/tex/works/Papers/2008/Pile_in_liquefied_soil_upU/SFSIModelSetup_upU_01.eps}}
{\includegraphics[width=9cm]{NewFiga/SFSIModelSetup_upU_01.eps}}
\caption{\label{coupling} Schematic description of coupling of displacement DOFs
($u_i$) of beam-column element (pile) with displacement DOFs ($u_i$) of u-p-U elements (soil).} 
\end{center}
\end{figure}
%



Nonlinear  force  based  beam--column  elements \citep{Spacone1996a, Spacone1996b} were used for modeling the 
pile--column. 
Pile  was assumed to be made of aluminum. This was done
in  order  to  be able to validate simulations with centrifuge experiments (when
they become available). Presented models were all done in prototype scale, while
for  (possible  future)  validation, select results will be carefully scaled and
compared    with   appropriate centrifuge modeling. Pile and the column were assumed to
have  a  diameter  of  $d=1.0$~m,    with Young's modulus of $E=68.5$~GPa, yield
strength  $f_y  =  255$~kPa,  and  the  density  $\rho  =2.7~{\rm kg/m^3}$. Wall
thickness   of  prototype  pile--column  is  $t  =  0.05$~m.  Lumped mass of  pile and
column    was distributed along the beam--column nodes, while an additional mass
was  added  on  top  $(m  =  1200$~kg)  that  represents  (small)  part  of  the
superstructure  mass.    This  particular  mass  $(m   =  1200$~kg) comes from a
standard (scaled up in our case) centrifuge model for pile--column--mass used at
UCD.

% 
%  %The critical state characteristics of the material model 
%  %and thin thickness soil elements around the pile hole can numerically 
%  %simulate the effects of gapping and slipary between the pile and soil \citep{Desai1984}.  
%  %The single pile is assumed made of aluminum with 
%  %the Young's modulus of 68.5 GPa,
%  %yield strength 255 kPa, 
%  %and the density 2.7 kg/$m^3$. 
%  %The pile cross section is an aluminum tub wrapped with a shrinkage layer. 
%  %The pile hole diameter is 1.2 m. 
%  %The aluminum tubbing has an outer diameter of 1.0 m and the wall thickness is 0.05 m. 
%  %In calculation, the aluminum tubing dimensions are used to define the pile fiber section. 
%  %The outer shrinkage layer is 0.1 m. 
%  %The mass of pile are distributed over nodes as lumped masses. 
%  
%  
%  % 
%  % In addition, at the top of the pile has a total structure mass of 1200 Mg. 
%  % The fiber section is defined by the OpenSees Steel01 material. 
%  % The pile section is subdivided into 20 segments in the circumferential direction and 10 segments in the radial direction. 
%  % The self-weight analysis of the pile is followed by soil self-weight analysis. 
%  % This analysis altogether with the previous stage 
%  % provide stable stress and pore pressure fields 
%  % with the pile inside the soils before the earthquake shaking.  
%  
Figure~(\ref{FEmesh3D}) (right side) shows side view of the column-pile-soil model
after second stage of loading.
%
%   \begin{figure}[!htbp]
%   \begin{center}
%   \includegraphics[width=0.6\textwidth]{NewFiga/FEmeshPileBeam.eps}
%   \caption{
%     \label{FEmeshPileBeam}  
%     {    Side view of the pile-soil model with  some element and node
%     annotation, used to visualize results.}
%   }  
%   \end{center}
%   \end{figure}
%   


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Third Loading Stage: Seismic Shaking}

After    the  application of self weight on the uniform soil profile, excavation
and  construction of the single pile with column and super structure mass on top
and  application  of  their self weight, the model is at the appropriate initial
state  for further application of loading. In this case, this additional loading
comprises  seismic  shaking.
%
%
%
For this stage, fixed horizontal DOFs used on the side planes
during the first stage are removed (set free).


The  input  acceleration  time  history, shown in Figure~(\ref{SinglePileGM}) 
was taken  from the recorded horizontal
acceleration  of Model No.1 of VELACS project \cite{Arulanandan93} by Rensselaer
Polytechnic  Institute  (\url{http::/geoinfo.usc.edu/gees/velacs/}).  
%
%%%%%
\begin{figure}[!htbp]
\begin{center}
\includegraphics[width=0.65\textwidth]{NewFiga/GM.eps}
\caption{
  \label{SinglePileGM}  
  {Input earthquake ground motions.}
}  
\end{center}
\end{figure}
%%%%%
%
The  magnitude  of the  motion is close to $0.2$~g, while   main shaking lasts
for about 12 seconds (from $1$~s to $13$~s).
%
Although   the  input  earthquake  motions  lasts  until  approx.  $13$~seconds,
simulations   are  continued  until  $120$~seconds  so  that  both  liquefaction
(dynamic)  and  pore  water  dissipation  (slow  transient) can be appropriately
simulated during and after earthquake shaking \citep{Jeremic2007e}.



% 
%  In addition to the horizontal level soil ground, 
%  cases with a small slope soil ground are investigated as well.  
%  This small slope plays an importance role during and after the earthquake excitation.
%  The effects of the small slope are approximately simulated 
%  by a horizontal distributing load loaded at the left side surface of the soil block, 
%  which is assumed 3\% of the soil vertical effective pressure at the same levels. 
%  This side load is added after the analysis of the pile driving and takes 
%  effects during and after the earthquake shaking.
%  The lateral load pushes the soils and causes the soil lateral spread.  
%  The soil lateral spread thus further pushes the pile and causes the pile deflects. 
%  The movement of the pile further results in the displacements of the upper structure. 
%  However, the real mechanism is not so simple during the earthquake,  
%  because there are interactions between the soils, pile and upper structure. 
%  it is also interesting to see the soil behaviors subjected to the earthquake loading
%  with different pile head conditions, namely, free head and fix head.  
%  The free pile head condition assumes the pile head have free horizontal and rotational DOFs 
%  in the lateral direction, which is to simulate condtionas 
%  similar to the bridge pier movement in the transverse direction.
%  The fixed pile head condition assumes the pile head have free horizontal but no rotational DOFs 
%  in the lateral direction, which is to mimic conditionas 
%  such as the bridge pier shaking in the longitude direction.  






%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Free Field, Lateral and Longitudinal Models}
\label{six_models}

Six      models    were  developed during the course of this study. First three models
(model numbers I,
II and III) were for level ground, while last three models (model numbers IV, V,
and  VI)  were  for  sloping ground. First in each series of models (model I for
level  ground  and  model  IV  for  sloping ground) were left without the second
loading  stage,  without  a  pile--column system. Other four models (numbers II,
III, V and VI)
were  analyzed  for  all  three  loading stages. Second in each series of models
(models  number  II  and  V)  had  all displacements and rotations of pile--column top 
(where additional mass representing superstructure was placed) left free,
without restraints. Thus, these two models represent lateral behavior of a bridge.
Third  in  each series of models (model numbers III and VI) had rotations in $y$
directions    fixed    at    the    pile--column top,  thus representing
longitudinal behavior of a bridge. Modeling longitudinal behavior of a bridge by
restraining  rotations perpendicular to the bridge superstructure is appropriate
if  the stiffness of a bridge superstructure is large enough, which in this case
it  was, as  it  was assumed to be a post--tensioned concrete box girder, so that realistically,
the top of a column does not rotate (much) during application of loads.
%
Table \ref{CasesOfSiglePile} summarizes models described above.
%
%
%%%%%
\begin{table}[!htbp]
\caption{ \label{CasesOfSiglePile}      
        Cases descriptions.  
}
\begin{center}
\begin{tabular}{|c|c|c|} 
\hline
      Case &  Model sketch & Descriptions     \\
\hline \hline
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      I    
      &
\includegraphics[width=0.025\textwidth,angle=0]{Model_I.eps}
      &   
      horizontal ground, no pile    \\ \hline 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      II   
      &
\includegraphics[width=0.025\textwidth,angle=0]{Model_II.eps}
      &   
      horizontal ground, single pile, free column head  \\ \hline 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      III  
      &
\includegraphics[width=0.025\textwidth,angle=0]{Model_III.eps}
      &   
      horizontal ground, single pile, no rotation at column head   \\
\hline
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      IV   
      &
\includegraphics[width=0.025\textwidth,angle=0]{Model_IV.eps}
      &   
      sloping ground, no pile      \\ \hline 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      V    
      &
\includegraphics[width=0.025\textwidth,angle=0]{Model_V.eps}
      &   
     sloping ground, single pile, free column head    \\ \hline 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      VI   
      &
\includegraphics[width=0.025\textwidth,angle=0]{Model_VI.eps}
      &   
      sloping ground, single pile, no rotation at column head    \\
\hline
\end{tabular}
\end{center}
\end{table}
%%%%%

% 
%  Totally six cases are listed in Table \ref{CasesOfSiglePile}. 
%  For ground conditions, 
%  cases of horizontal ground without any slope and small slope ground with a slope angle 3\% are considered. 
%  For pile head conditions, 
%  cases of no pile, single pile with free pile top, and single pile without rotation at the top are taken into account. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Simulation Results}
\label{Single_Pile_Liquefiable_Soil}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\subsection{Pore Fluid Migration}
\label{Numerical_Results} 


Figures  (\ref{Ru_14}) through (\ref{Ru_36}) show the Ru time history for  up to
$30$~seconds, for elements (at one of  Gauss point) e1, e3, e5 and e7 (refer to
right side of Figure~(\ref{FEmesh3D})).
%
It    is    important  to  note  that $R_u$ is defined as the ratio of the
difference  of initial mean and current mean effective stresses over the initial mean
effective stress:
%
\begin{equation}
R_u = \frac{p_{initial}^{\prime} - p_{current}^{\prime}}{p_{initial}^{\prime}}
\nonumber
\end{equation}
%
where  mean effective
stress is defined as $p^{\prime} = \sigma^{\prime}_{kk}/3$.
%
This  is  different  from traditional  definition for $R_u$, that uses ratio of 
excess pore pressure
over the initial mean
effective   stress   ($p_{initial}^{\prime}$).   However, these two  definitions  are
essentially    equivalent,  as soil is  in  the state of liquefaction for
$R_u=1$ (so that $p_{current}^{\prime}=0$), 
while there is no excess pore pressure for $R_u=0$ (so that 
$p_{initial}^{\prime} = p_{current}^{\prime}$).
%
However, the former
definition  is  advocated  here as it  avoids the interpolation of 
pore    pressure  or  extrapolation  of  the  stresses (as the latter definition
requires), since for  the  u-p-U  element, stresses are available at Gauss
points while pore pressures are available element nodes.
%
\begin{figure}[!htbp]
\begin{center}
\includegraphics[width=0.8\textwidth]{NewFiga/Ru_case_1_4.eps}
\vspace*{-0.6cm}
\caption{\label{Ru_14} $R_u$ times histories for elements e1 (top element), 
e3, e5, and e7 (bottom element)
Gauss  point)  for  Cases  I  (level ground, no pile) and IV (sloping ground, no
pile).}
\end{center}
\end{figure}
%
In  particular,  Figure  \ref{Ru_14} shows $R_u$ time histories for four points for
models I (level ground without pile) and model IV (sloping ground without pile).
%
It   is  noted  that  differences  are  fairly  small. 
%
It is interesting to observe  that  lower  layers  do not liquefy as supply
of  pore  fluid  for  initial void ratio of $e_0=0.8$ is too small, and the pore
fluid dissipation upward seems to be to rapid.
%
On the other hand, the upper soil layers do reach close to or
liquefaction  state  ($R_u=1$).  
%
This is primarily due to the propagation of pore fluid  pressure/volume  from   
lower layers upward (pumping effect) and, in addition to that, to a local excess
pore fluid production.
%
These results can also be
contrasted  with  those  of \citet{Jeremic2007e}, where similar pumping scenario
has been observed. 
%
The main difference between soil
used  by  \citet{Jeremic2007e}  and  here is in the coefficient of permeability.
%
Namely,      here     $k     =     1.0\times     10^{-4}$~$m/s$     was     used
\citep{Cubrinovski2008,Uzuoka2008}         while
\citet{Jeremic2007e} used $k = 5.0\times 10^{-4}$~$m/s$. 
%
It  is  important  to note that   other  values  of  permeability  for  Toyoura 
sand  have  also  been reported \citep{Sakemi1995}, but current value was chosen
based on \citet{CubrinovskiPrivateCommunications}.
%

In  addition  to  that, similar to \citet{Jeremic2007e}, sloping ground case shows
larger  $R_u$  spikes,  since  there is static shear force (stress) that is always
present from gravity load on a slope. This static gravity on a slope will result 
in an asymmetric
horizontal shear stresses in the down--slope direction during cycles of
shaking. 
%
This asymmetric shear stress induces a more dilative response for down--slope 
shaking which will help soil regain its
stiffness  in  the dilative parts of the loading cycles. This observation can be
used to explain smaller $R_u$ spikes for the sloping ground case.
%
Of  course,  this  asymmetry  in  loading  will  result  in larger accumulation 
of down--slope deformation.
%




While  $R_u$  ratios  for  level  and  sloping ground cases are fairly similar
along  the  depth  of  the model, the response changes when the pile is present.
Figure~(\ref{Ru_25}) shows $R_u$ responses at four different points (along the depth)
approximately  midway  between  the pile and the model boundary, in the plane of
shaking  (see  location  of  those  elements  in Figure~(\ref{FEmesh3D}) on page
\pageref{FEmesh3D}).
%
In  comparison  to  behavior  without  the  pile  (Figure~(\ref{Ru_14})), 
 it  is  immediately  obvious  that  addition of a pile 
with a mass on top reduces
$R_u$  during  shaking for the top element (e1). This is to be expected as presence
of  a  pile--column--mass (PCM) system changes the dynamics of the top layers of soil
significantly  enough to reduce total amount of shear. This is particularly true
for    the  top  layers  of  soil  as  effects  of  column--mass  tend to create
compressive  and  extensive  movements  (compression when the PCM
system  moves  toward  soil  and  extension, and possibly even tension, when PCM
system    moves    away   from  soil).  However, this  extension,  or  possible 
tension, is not
directly  observable  in presented plots since array of elements where we follow
$R_u$  (e1,  e3,  e5,  e7) is some distance away from the pile--soil interface.
% and possibly forming
%tension gap between pile and adjacent soil elements.
%
Middle layers (e3 and e5), on the other hand, display similar response to that of
Cases  I and IV, as shown in Figure~(\ref{Ru_14}). It is
noted  that  in a case with of sloping ground with pile, the $R_u$ measurements are
always  larger  that those for level ground (this is also observed for Cases III
and  VI,  as  shown  in  Figure~(\ref{Ru_36})). This is
expected as presence of a pile in loose sand, and particularly the dynamic movement
of  a  pile  during  seismic  shaking, create an additional shearing deformation
field  (in  the  soil  adjacent  to  the  pile)  that  provides  for  additional
(faster) compression  of  soil  skeleton  and  thus creates additional volume of pore
fluid, that is then distributed to adjacent soil  (adjacent to the pile).
%
%
\begin{figure}[!htbp]
\begin{center}
\includegraphics[width=0.8\textwidth]{NewFiga/Ru_case_2_5.eps}
\vspace*{-0.6cm}
\caption{\label{Ru_25} $R_u$ times histories for elements e1, e3, e5, and e7 (upper
Gauss  point)  for  Cases  II  (level  ground, with pile--column, free column head) and V
(sloping ground, with pile--column, free column head).}
\end{center}
\end{figure}
%

Particularly  interesting  are $R_u$ results  for soil element e7, which is
located   below   pile   tip   level   (see  Figure~(\ref{FEmesh3D})). 
Observed  $R_u$  for Case V in element e7 is significantly larger than for
the same element for Case  II.  Similarly,  simulated $R_u$ is larger than what was observed 
in cases without a pile  (see  bottom  of  Figure~(\ref{Ru_14})).  This
increase  in   $R_u$  for  Case  V  (sloping ground with pile) is
explained    by    noting    that  the pile ``reinforces" upper soil layers 
and thus  prevents
excess    shear deformation in the upper $12.0$~m  of soil (above pile  tip).
% ,  when  compared with sloping ground Case without the pile (IV).
%This  reduced  deformation  (and subsequently shearing). 
The reduction of  deformation
in  upper  layers of soil (top $12.0$~meters) results in transfer of excessive soil deformation
to        soil layers below pile tip (where element e7 is located). This, in
turn, results  in a much  larger and faster shearing of those lower loose soil layers. This
significantly    larger   shearing results in a much higher $R_u$. Deformed shape,
shown    in    Figure~(\ref{Snapshots123456})  for  Case  V,  reinforces  this
explanation, showing much
large   shearing  deformation  in  lower  soil  layers,  below  pile  tip.  Same
observation        can    be    made    for    Case    VI,    shown    in  
Figure~(\ref{Snapshots123456}).
%


Observation  similar  to the above, for Cases II and V can be made for Cases III
and VI, results for which are shown in Figure~(\ref{Ru_36}).
%
\begin{figure}[!htbp]
\begin{center}
\includegraphics[width=0.8\textwidth]{NewFiga/Ru_case_3_6.eps}
\vspace*{-0.6cm}
\caption{\label{Ru_36} 
$R_u$ times histories for elements e1, e3, e5, and e7 (upper
Gauss  point)  for  Cases  III  (level  ground, with pile, no rotation of pile
head) and VI (sloping ground, with pile-column, no rotation of column head).}
\end{center}
\end{figure}
%
One       noticeable  difference in $R_u$ results between cases with free column
head (Cases II and V) and  cases with fixed rotation
column head (Cases III and VI) is in    significantly   
larger  (and  faster) development  of $R_u$ close to soil surface for a stiffer, no
rotation  column  cases (Cases III and VI).
% Similar observation is made for sloping ground
%cases (Cases V, sloping ground with pile, restricted rotation on
%top     of  column) and VI (sloping ground with pile, restricted rotation on top
%of column). 
%
This  much  larger  $R_u$  observed in a  ``stiffer" PCM system setup, is
due to larger shearing deformation that develops in soils adjacent to the
pile    during    shaking.  The stiffer PCM system can displace less (because of
additional no rotation condition on column top) while the
soil   beneath  is undergoing shaking (same demand in all cases), thus resulting
in  larger  relative  shearing  of soil, which then results in larger and faster
pore  pressure  development  close  to  the soil surface, where the column no
rotation effect is most pronounced.






%   %  
%   In order to compare the slope effects on $R_u$, 
%   $R_u$ for 
%   Case I and IV (Figure \ref{Ru_14}), 
%   Case II and V (Figure \ref{Ru_25}), 
%   Case III and VI (Figure \ref{Ru_36})are plotted together up to 30 seconds.. 
%   Form these figures, $R_u$ at all elements are found to increase shortly after earthquake shaking. 
%   At upper layers, $R_u$ is up to as big as close to 1.0, 
%   which means that liquefaction have occurred at upper layers. 
%   The $R_u$ time histories are almost overlapped for Case I and IV, 
%   which indicates that the slope effect on $R_u$ is not distinct 
%   for free field without pile (Figure \ref{Ru_14}). 
%   However, $R_u$ curves with level ground are generally lower than those with slope ground, 
%   which suggests strong slope effect on $R_u$ for the cases with the pile in soil, 
%   whichever the pile head is fixed or not (Figure \ref{Ru_25} and Figure \ref{Ru_36}).
%    


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

\subsection{Soil Skeleton Deformation}
%%


A  number  of  deformation  modes  is observed for both level and sloping ground
cases,   with   or  without  PCM  system.  Figure~(\ref{Snapshots123456})  shows
deformation  patterns  and  excess  pore pressures in symmetry plane for all six
cases over a period of eighty seconds.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[!htbp] 
\begin{center}
\begin{tabular}{rlllllll} 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
I
&
\includegraphics[width=0.04\textwidth,angle=0]{Model_I.eps}
&
\includegraphics[height=0.12\textwidth]{NewFiga/Snap_1_T002.eps}
&
\includegraphics[height=0.12\textwidth]{NewFiga/Snap_1_T005.eps}
&
\includegraphics[height=0.12\textwidth]{NewFiga/Snap_1_T010.eps}
&
\includegraphics[height=0.12\textwidth]{NewFiga/Snap_1_T015.eps}
&
\includegraphics[height=0.12\textwidth]{NewFiga/Snap_1_T020.eps}
&
\includegraphics[height=0.12\textwidth]{NewFiga/Snap_1_T080.eps}
\\
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
II
&
\includegraphics[width=0.04\textwidth,angle=0]{Model_II.eps}
&
\includegraphics[height=0.12\textwidth]{NewFiga/Snap_2_T002.eps}
&
\includegraphics[height=0.12\textwidth]{NewFiga/Snap_2_T005.eps}
&
\includegraphics[height=0.12\textwidth]{NewFiga/Snap_2_T010.eps}
&
\includegraphics[height=0.12\textwidth]{NewFiga/Snap_2_T015.eps}
&
\includegraphics[height=0.12\textwidth]{NewFiga/Snap_2_T020.eps}
&
\includegraphics[height=0.12\textwidth]{NewFiga/Snap_2_T080.eps}
\\
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
III
&
\includegraphics[width=0.04\textwidth,angle=0]{Model_III.eps}
&
\includegraphics[height=0.12\textwidth]{NewFiga/Snap_3_T002.eps}
&
\includegraphics[height=0.12\textwidth]{NewFiga/Snap_3_T005.eps}
&
\includegraphics[height=0.12\textwidth]{NewFiga/Snap_3_T010.eps}
&
\includegraphics[height=0.12\textwidth]{NewFiga/Snap_3_T015.eps}
&
\includegraphics[height=0.12\textwidth]{NewFiga/Snap_3_T020.eps}
&
\includegraphics[height=0.12\textwidth]{NewFiga/Snap_3_T080.eps}
\\
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
IV
&
\includegraphics[width=0.04\textwidth,angle=0]{Model_IV.eps}
&
\includegraphics[height=0.12\textwidth]{NewFiga/Snap_4_T002.eps}
&
\includegraphics[height=0.12\textwidth]{NewFiga/Snap_4_T005.eps}
&
\includegraphics[height=0.12\textwidth]{NewFiga/Snap_4_T010.eps}
&
\includegraphics[height=0.12\textwidth]{NewFiga/Snap_4_T015.eps}
&
\includegraphics[height=0.12\textwidth]{NewFiga/Snap_4_T020.eps}
&
\includegraphics[height=0.12\textwidth]{NewFiga/Snap_4_T080.eps}
\\
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
V
&
\includegraphics[width=0.04\textwidth,angle=0]{Model_V.eps}
&
\includegraphics[height=0.12\textwidth]{NewFiga/Snap_5_T002.eps}
&
\includegraphics[height=0.12\textwidth]{NewFiga/Snap_5_T005.eps}
&
\includegraphics[height=0.12\textwidth]{NewFiga/Snap_5_T010.eps}
&
\includegraphics[height=0.12\textwidth]{NewFiga/Snap_5_T015.eps}
&
\includegraphics[height=0.12\textwidth]{NewFiga/Snap_5_T020.eps}
&
\includegraphics[height=0.12\textwidth]{NewFiga/Snap_5_T080.eps}
\\
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
VI
&
\includegraphics[width=0.04\textwidth,angle=0]{Model_VI.eps}
&
\includegraphics[height=0.12\textwidth]{NewFiga/Snap_6_T002.eps}
&
\includegraphics[height=0.12\textwidth]{NewFiga/Snap_6_T005.eps}
&
\includegraphics[height=0.12\textwidth]{NewFiga/Snap_6_T010.eps}
&
\includegraphics[height=0.12\textwidth]{NewFiga/Snap_6_T015.eps}
&
\includegraphics[height=0.12\textwidth]{NewFiga/Snap_6_T020.eps}
&
\includegraphics[height=0.12\textwidth]{NewFiga/Snap_6_T080.eps}
\\
&
t=
&
2~sec
&
5~sec 
&
10~sec
&
15~sec
&
20~sec
&
80~sec
\end{tabular}
\includegraphics[width=10cm,height=1cm]{NewFiga/GMklot02.ps}\hspace*{1cm}
\\
\includegraphics[angle=-90,width=0.5\textwidth]{NewFiga/Snap_scale.ps}
\caption{
  \label{Snapshots123456}
  {Time  sequence  of  deformed  shapes and excess  pore  pressure in symmetry  
  plane  of  a soil system. Deformation is exaggerated 15 times; Color scale for
  excess    pore    pressures  (above)  is  in $kN/m^2$. Graph of ground motions
  used   (also   shown   in  Figure~(\ref{SinglePileGM})) is placed  below
  appropriate  time snapshots and is matching for $t=2, 5, 10, 15, 20$~seconds while
  at $t=80$~seconds there is no seismic shaking.}
}
\end{center}
\end{figure}
%
A    number    of  observation can be made on both deformation patterns,  excess
pore fluid patterns and their close coupling.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\paragraph{Level   Ground     without  Pile (Case I).}  
%
Excess  pore  pressures  and  deformations in symmetry plane  for  level   ground   without    a  
pile are shown in Figure~(\ref{Snapshots123456}) (I).
%
At  the very beginning (at $t=2$~s) there is initial
development  of excess pore fluid pressure in the middle soil layers. This   
expected,   as  the self weight loading stage has densified lower soil
layers enough so that their response is not initially contractive enough  to produce
excess pore pressure. Top soil layers, on the other hand, have a drainage boundary
(top  surface)  too  close  to  develop any significant excess pore pressures. As
seismic  shaking  progresses  (for  $t  =  5,  10$~s),  the  excess pore pressure
increases,  and  starts  developing  in  lower soil layers as well. It should be
noted that a  small non-uniformity  in  results  is present. For example, zones of
variable, nonuniform excess pore pressures on the lower mid and right side for Case I at
$t=10$~s  develop.    
%This  non-uniformity was not initially expected, as the
%boundary  conditions  and  the model should be similar to a 1D soil column (Case I).
%However,  nonuniform mesh (many small, long elements in the middle,
%large  elements  outside  this  middle  zone)  introduce small errors in results
%(which  can be observed by slightly nonuniform results at $t=12$~s and $t=10$~s).
%%
Nonuniform mesh (many small, long elements in the middle,
large  elements  outside  this  middle  zone) may introduce small numerical errors in results
which  can be observed by slightly nonuniform results at $t=10$~s and $t=15$~s.
%%
It  should  be  noted  that results for excess pore pressure shown for
first  $13$~seconds  (during shaking) in Figure~(\ref{Snapshots123456}) (I) are
transient in nature, that is,
seismic  waves are traveling throughout the domain (model) during shaking (first
$13$~seconds)  and  slight  oscillations  in  vertical stresses are possible. This
oscillations will contribute to the (small) non-uniformity of excess pore pressure 
results.  After  the shaking (after $15$~seconds) resulting excess pore pressure
field is quite uniform.



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\paragraph{Level   Ground     with  Pile (Cases II and III).}  
%
Excess  pore  pressures  and  deformations in symmetry plane for models with PCM
system  and with two different boundary conditions at
top  of  column (see  model description in section~\ref{six_models}) 
in  level   ground  are shown   in       Figures~(\ref{Snapshots123456})  (II and III).  
%
One of the interesting observations is  significant shearing
and    excess  pore pressure generation adjacent to the pile tip. The reason for
this is that   pile  is  too  short,  that  is,  pile tip has significant
horizontal    displacements  during  shaking.  Those  pile  tip  displacements  
shear the soil, resulting in excess pore pressure generation.
As  soon  as  there is time for dissipation, this localized excess pore pressure
dissipates to adjacent soil, and then, after shaking has ceased (at $t=13$~s and
later), it slowly dissipates upward.
%
Addition  of  pile  into  the  model (construction), with a highly impermeable
elements  (that  mimic  permeability  of concrete) is apparent as there is a low
 excess pore pressure region in the middle of model, where pile is located.




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\paragraph{Sloping  Ground     without  Pile (Case IV).}  
%
Excess  pore  pressures  and  deformation  in  symmetry plane for sloping ground
without pile is shown in Figures~(\ref{Snapshots123456}) (IV).
%
It  is  noted  that  initially the excess pore pressure starts developing in 
middle    soil    layers,similar  to  the  Case  I above.  Bottom  layers  start
developing excess pore
pressure  only  after  significant shear deformation occurs (at
$t=10$~s)  at approximately $2/3$ of the model depth. Lower layers have
densified enough during self weight stage of loading that initial shaking is not
strong enough to create excess pore water pressure, rather, those layers are fed
by  the excess pore pressure from above. Lower soil layers also do not develop much
deformation,   while  middle  and  upper  layers    together  develop  excessive
horizontal deformation.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\paragraph{Sloping  Ground  with  Pile (Cases V and VI).}  
%
Excess  pore  pressures  and  deformation  in  symmetry plane for sloping ground
with PCM system are shown in Figures~(\ref{Snapshots123456}) (V and VI).
%
Similar to the above cases (II and III), pile is too short  and there 
is  again  excessive shearing of soil at the pile tip, suggesting large movement
of that pile tip.
%
In  addition to that, pile introduces significant stiffness to upper $12$~meters
of  soil  (along the length of pile) and helps reduce deformation of those upper
soil  layers.  Down-slope  gravity  load  is thus  transferred to lower soil layers
(below pile tip) which exhibit most of the deformation. It should be noted that
soil  in  middle and upper  layers  (adjacent to pile) does deform, just not as much as the
soil below pile tip. The predominant mode of deformation of middle  soil layers is
shearing    in  horizontal  plane,  around  the  pile. Deformation in horizontal
plane  is  not  significant  as the pile is short in this examples (as mentioned
above)  and  does  not have enough horizontal support at the bottom. 
The  deformation pattern of a  soil -- pile system is such that pile experiences
significant rotation, and deforms with the soil that moves down-slope.
%
If  the  pile  was  longer,  and if it had significant horizontal support at the
bottom, the middle and upper soil layers would have showed more significant flow
around the pile in horizontal planes.

%
Upper     layers     undergo     significant     settlement,    as    seen    in
Figure~(\ref{Settlements}).
This  settlement  is mainly caused by the above mentioned rotation of pile--soil
system, where soil in general settles (compacts) but also undergoes differential
settlement, between left (up--slope from pile) and right (down--slope from
pile) side of the model.
%
As significant  shearing  with excess pore pressure generation develops in lower soil
layers,  below  pile  tip,  those lower layers contribute to most of down--slope
horizontal deformation.
%
In a sense, all the demand from down-slope gravity forces
and  seismic  shaking  is  now  responded to by lower  soil layers,
which   contribute   to   most  of  the  excess  pore  pressure  generation  and
consequently, to most of the soil deformation. 
%
Soil  surface  horizontal deformation is thus strongly influenced by significant
shearing  of  the  bottom  layers  and  by rotation of the middle and upper soil
layers with the pile.
%
%
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[!htbp]
\begin{center}
\begin{tabular}{rrr} 
%\hline 
I
&
II
&
III
\\
\includegraphics[height=0.12\textwidth]{NewFiga/settlement_case1.ps}
&
\includegraphics[height=0.12\textwidth]{NewFiga/settlement_case2.ps}
&
\includegraphics[height=0.12\textwidth]{NewFiga/settlement_case3.ps} 
\\
\\
% \hline
IV
&
V
&
VI
\\
\includegraphics[height=0.12\textwidth]{NewFiga/settlement_case4.ps}
&
\includegraphics[height=0.12\textwidth]{NewFiga/settlement_case5.ps}
&
\includegraphics[height=0.12\textwidth]{NewFiga/settlement_case6.ps} 
%\\
%
%\hline
\end{tabular}
\includegraphics[angle=-90,width=0.50\textwidth]{NewFiga/settlement_scale.ps}
\caption{
  \label{Settlements}
{Soil  surface    settlements at $120$~s for all six cases. Color scale given in
meters}
}
\end{center}
\end{figure}
%
%
It  is  interesting  to  note  that  the  largest  settlement  is  observed just
down-slope from pile for Cases V and VI.

% 
%  It  is  also  very  interesting  to  note  that, even though soil achieves cyclic
%  mobility      (for    example  for  Cases  II  and  III,  see  $R_u$  values  in
%  Figures~(\ref{Ru_25}), 
%  and  (\ref{Ru_36}))  and  in  some  cases  full
%  liquefaction   (for   example  for  Cases  V  and  VI,  see  $R_u$  values    in
%  Figures~(\ref{Ru_25})
%  and   (\ref{Ru_36})), pile is still able to
%  hold  the  vertical  loads (from mass on top and from the its own and
%  column's self weights) by utilizing skin friction.
%  

% 



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

\subsection{Pile Response}



%
%
% 
%   for plots the bending moment distribution at each calculation
%  time step for Case II, III, V and VI.
%  During the structure-pile-foundation interaction, 
%  the moment envelope upon the time duration is important to pile performance. 
%  The pile moment distributions between the free and fixed pile head are quite different.  
%  For free pile head, the maximum pile bending moments are at depth between 0 and -2 feet 
%  with a magnitude about 800 kN, 
%  while for fixed pile head, the maximum pile bending moment are at the pile head 
%  with a magnitude close to 2000 kN. 
%  The slope effect on the bending moment envelop is not distinct. 
%  

Figure~(\ref{Moments})    shows  bending  moment envelops for pile--column--mass
(PCM) system for all four cases (II, III, V and VI). 
%
It  should  be noted that bending moment diagrams are plotted on compression side
of the beam--column.
%
A number of observations can be made about bending moment envelopes.
%
For  cases  with  free pile head (shaking transverse to the bridge main axes,
Cases II and V)
the  maximum  moments  are  attained  in soil, at depths of approximately $0.6D -
1.2D$,    where    $D$~($=1.0$~m in this case)  is  the pile diameter.
%
Opposed to that are cases for PCM systems with restricted rotations at
the   pile top which (Cases III and VI), which, of course feature largest moment
at the column top.
Maximum   bending  moments for section of PCM system in soil (pile) in these two
cases are now attained at the depth of approximately $1.8D - 2.0D$.



It   is  noted  that  bending  moment  envelopes  are  mostly  symmetric.  Slight
non--symmetry is introduced for cases on sloping ground (Case V and VI).
It  is  also  noted that moments do exist (are not zero) all the say to the bottom of the pile.
Theoretically, moments should be zero at the pile tip, but since physical volume
of the pile is considered (see note on that in section \ref{pile_installation}
and  Figure~(\ref{coupling})),   differential pressure on pile bottom from soil
will  produce  small (non--zero) moments even at the pile tip. More importantly,
non--zero  moments  at the bottom and along the lower part of the pile show that
pile  is  indeed  too  short,  and thus changing curvatures are present along the
whole length of the pile.


%
\begin{figure}[!htbp]
%\vspace*{-1cm}
\begin{center}
%\begin{tabular}{rrrr} 
\begin{tabular}{rr} 
%II
%&
\includegraphics[width=0.3\textwidth]{NewFiga/M/M2.eps}
&
%III
%&
\includegraphics[width=0.3\textwidth]{NewFiga/M/M5.eps}
\\
%V
%&
\includegraphics[width=0.3\textwidth]{NewFiga/M/M3.eps}
&
%VI
%&
\includegraphics[width=0.3\textwidth]{NewFiga/M/M6.eps}
\end{tabular} 
\caption{
  \label{Moments}
{Envelope   of    bending  moments  for pile--column system for Cases II, III, V
and VI.}
}
\end{center}
\end{figure}



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

\subsection{Pile Pinning Effects}


Piles  in sloping liquefying ground can also be used to resists movement of soil
(all  liquefied  or  liquefied  with  hard crust on top) down--slope. For models
developed  in  this paper, pile pinning effect can be investigated for
Cases  IV,  V  and  VI. In particular, deformation of sloping ground without the
pile  (Case  IV)  can  be  compared with either of the cases of piles in sloping
ground, Cases V and VI.
%
It  is  very important to note, again, that models developed here had relatively
short  pile,  and that major soil shearing developed below the pile
tip.  This  apparent  shortcoming  of a short pile results in reduced pile pinning
capacity,    thus   reducing  the  down--slope  movement  by only approximately 
half, from
$0.35$~m (Case IV)   to    $0.22$~m   (Case  V)  and to  $0.18$~m (Case VI) as seen
in Figure~(\ref{Dx}).  It would have
been  expected  that, had the pile been longer and had it penetrated in deeper,
non-liquefiable layers, it would have reduced down--slope movement of the soil to
a    much  larger  extent.  However,  had  the  pile  been longer and had it penetrated
non-liquefiable layers, it would have had a much firmer horizontal support at the
bottom and would have thus attracted much larger forces too, potentially leading
to pile damage and yielding.
%  % 
%  OVDE 
%  
%  
%  Figure \ref{Dx} plots the lateral displacements of the ground center
%  (the middle point at the ground surface in the middle plain, 
%  or the point of pile at the ground surface). 
%  For all cases, the ultimate lateral displacements of the ground center, 
%  there are some residual lateral displacement after the earthquake exciting. 
%  It is obviously that the pile has much bigger lateral displacement for the slope ground (Case I, II and III)
%  than those for the level ground (Case IV, V, and VI), which is as expected. 
%  No matter for level ground, or for slope ground, 
%  the the lateral displacements of the ground center without pile (Case I and IV) are much higher 
%  than those with pile (Case II, III and V and VI), 
%  which suggest strong pile pinning effects on the ground lateral displacement. 
%     
%  
\begin{figure}[!htbp]
\begin{center}
%\includegraphics[width=0.6\textwidth]{NewFiga/Dx123.eps}
%\\
\vspace*{1cm}
\includegraphics[width=0.7\textwidth]{NewFiga/Dx456.eps}
\\
\caption{
  \label{Dx}
  {Down--slope   movement  at the ground surface (model center) for Cases IV (no
  pile), V and VI (with pile--mass system).}
}
\end{center}
\end{figure}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Summary}
\label{Summary}




Presented in this paper was methodology for numerical modeling and simulation of
piles  in  liquefiable soil. Of particular interest was the detailed description
of modeling which aimed at replicating the prototype model as close as possible.
High  fidelity  modeling included use of verified and validated models, detailed
model  development,  including  use  of realistic loading stages. 
%
Detailed  application  of  loading staged, starting from a  zero state of stress
and  strain  for  a  soil  without  a pile, followed by application of soil self
weight,   excavation   and   pile--column   installation   with  application  of
pile--column self-weighting is finally followed by seismic loading with extended
time after that for dissipation of excess pore pressures that have developed.
%
An  implementation  for  a  bounding  surface elastic--plastic  sand  model that 
accounts for fabric change,  and for a fully coupled porous media
(soil   skeleton)  -- pore fluid (water) dynamic finite element formulation were
developed and used in simulation  of soil and water displacement and pore
water  pressure.
%  


Six  models  were  developed  and  simulated,  feature  level and sloping ground
without and with pile--column systems.
%
Results of simulations are presented with the aim of increasing our 
understanding    of    behavior    of    soil--pile--column    systems    during
liquefaction  events,  including  lateral  soil  deformation,  effects  of  pile
pinning, and  ground settlement.
%
In  addition  to detailed presentation of useful and interesting results, one of
the  main  aims  of  this  paper  was  to emphasize the need for, importance and
availability of high fidelity modeling tools for simulating effects of liquefied
soil on soil--structure systems.




%occurring,  excess pore water building-up and dissipation, soil skeleton lateral
%deformation,  ground  settlement,  pile response, and pile-soil interaction were
%interpreted at some length.


  
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\noindent
\section*{Acknowledgment}
%{\Large \bf Acknowledgment}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{samepage}
%
Work  presented  in  this  paper  was  supported  in part by the Earthquake
Engineering  Research  Centers  Program of the National Science Foundation under
Award Number EEC-9701568.
%
\end{samepage}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\renewcommand{\baselinestretch}{1.1}
%\small\normalsize  % trick from Kopka and Daly p47.
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


\bibliography{refmech,refcomp,zhao,Mahdi,auxSFSI,/home/jeremic/tex/works/Proposals/2000/CalTrans/aux}
%\bibliography{refmech,refcomp,zhao,Mahdi,auxSFSI}
\bibliographystyle{unsrtnat}


\end{document}
\bye



