%% 
%% Copyright 2007, 2008, 2009 Elsevier Ltd
%% 
%% This file is part of the 'Elsarticle Bundle'.
%% ---------------------------------------------
%% 
%% It may be distributed under the conditions of the LaTeX Project Public
%% License, either version 1.2 of this license or (at your option) any
%% later version.  The latest version of this license is in
%%    http://www.latex-project.org/lppl.txt
%% and version 1.2 or later is part of all distributions of LaTeX
%% version 1999/12/01 or later.
%% 
%% The list of all files belonging to the 'Elsarticle Bundle' is
%% given in the file `manifest.txt'.
%% 

%% Template article for Elsevier's document class `elsarticle'
%% with numbered style bibliographic references
%% SP 2008/03/01
%%
%% 
%%
%% $Id: elsarticle.cls,v 1.20 2008-10-13 04:24:12 cvr Exp $
%%
%%
%\documentclass[preprint,12pt]{elsarticle}

%% Use the option review to obtain double line spacing
\documentclass[preprint,review,12pt]{elsarticle}

%% Use the options 1p,twocolumn; 3p; 3p,twocolumn; 5p; or 5p,twocolumn
%% for a journal layout:
%% \documentclass[final,1p,times]{elsarticle}
%% \documentclass[final,1p,times,twocolumn]{elsarticle}
%% \documentclass[final,3p,times]{elsarticle}
%% \documentclass[final,3p,times,twocolumn]{elsarticle}
%% \documentclass[final,5p,times]{elsarticle}
%% \documentclass[final,5p,times,twocolumn]{elsarticle}

%% if you use PostScript figures in your article
%% use the graphics package for simple commands
\usepackage{graphics}
%% or use the graphicx package for more complicated commands
\usepackage{graphicx}
%% or use the epsfig package if you prefer to use the old commands
\usepackage{epsfig}

%% The amssymb package provides various useful mathematical symbols
\usepackage{amssymb}
%% The amsthm package provides extended theorem environments
\usepackage{amsthm}

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

%% natbib.sty is loaded by default. However, natbib options can be
%% provided with \biboptions{...} command. Following options are
%% valid:

%%   round  -  round parentheses are used (default)
%%   square -  square brackets are used   [option]
%%   curly  -  curly braces are used      {option}
%%   angle  -  angle brackets are used    <option>
%%   semicolon  -  multiple citations separated by semi-colon 
%%   colon  - same as semicolon, an earlier confusion
%%   comma  -  separated by comma
%%   numbers-  selects numerical citations
%%   super  -  numerical citations as superscripts
%%   sort   -  sorts multiple citations according to order in ref. list
%%   sort&compress   -  like sort, but also compresses numerical citations
%%   compress - compresses without sorting
%%
%% \biboptions{comma,round}

% \biboptions{}


%Added by Mahdi
\usepackage{geometry}
\geometry{left=2.5 cm,right=2.5 cm,top=2.5 cm,bottom=2.5 cm}
\usepackage{amsfonts}
\usepackage{amsmath}
\def\SingleSpaced{\def\baselinestretch{1.1}\large\normalsize}
\def\DoubleSpaced{\def\baselinestretch{1.6}\large\normalsize}
\usepackage{subfigure}
\graphicspath{
{./Figures/}
}

\journal{Soil Dynamics and Earthquake Engineering}

\begin{document}

\begin{frontmatter}

%% Title, authors and addresses

%% use the tnoteref command within \title for footnotes;
%% use the tnotetext command for the associated footnote;
%% use the fnref command within \author or \address for footnotes;
%% use the fntext command for the associated footnote;
%% use the corref command within \author for corresponding author footnotes;
%% use the cortext command for the associated footnote;
%% use the ead command for the email address,
%% and the form \ead[url] for the home page:
%%
%% \title{Title\tnoteref{label1}}
%% \tnotetext[label1]{}
%% \author{Name\corref{cor1}\fnref{label2}}
%% \ead{email address}
%% \ead[url]{home page}
%% \fntext[label2]{}
%% \cortext[cor1]{}
%% \address{Address\fnref{label3}}
%% \fntext[label3]{}

\title{Propagation of Seismic Waves through Liquefied Soils}

%% use optional labels to link authors explicitly to addresses:
%% \author[label1,label2]{<author name>}
%% \address[label1]{<address>}
%% \address[label2]{<address>}

\author[1]{Mahdi Taiebat\corref{cor1}}
\ead{mtaiebat@civil.ubc.ca}
\author[2]{Boris Jeremi{\'c}}
\author[2,3]{Yannis F. Dafalias}
\author[4]{Amir M. Kaynia}
\author[5]{Zhao Cheng}

\cortext[cor1]{Corresponding author. Tel.: +1 (604) 822-3279; Fax: +1 (604) 822-6901.}

\address[1]{Department of Civil Engineering, University of British Columbia, Vancouver, B.C. V6T 1Z4, Canada}
\address[2]{Department of Civil and Environmental Engineering, University of California, Davis, CA 95616, USA}
\address[3]{Department of Mechanics, National Technical University of Athens, Zographou 15780, Hellas}
\address[4]{Norwegian Geotechnical Institute, P.O. Box 3930 Ullevaal Stadion, N-0806 Oslo, Norway}
\address[5]{Earth Mechanics Inc. Oakland CA 94621, USA}

\begin{abstract}
%% Text of abstract
For prediction of the earthquake response of saturated
porous media it is essential to correctly simulate the
generation, redistribution and dissipation of excess pore
water pressure during and after earthquake shaking.
%
To this end, a reliable numerical tool requires a dynamic,
fully coupled formulation for solid--fluid interaction and a
versatile constitutive model.
%
Presented in this paper is a three dimensional finite
element framework which has been developed and utilized for this
purpose.
%
The framework employs fully coupled dynamic field
equations with a $u-p-U$ formulation for simulation of pore
fluid and soil skeleton interaction, and a SANISAND
constitutive model for solid skeleton.
%
Following a detailed verification and validation of the
formulation and implementation of the developed numerical
tool, the numerical tool has been employed in the seismic response of saturated porous media. The study includes the mechanism of propagation of the earthquake-induced shear waves and liquefaction phenomena in uniform and layered profiles of saturated sand deposits.

\end{abstract}

\begin{keyword}
%% keywords here, in the form: keyword \sep keyword
%% PACS codes here, in the form: \PACS code \sep code
%% MSC codes here, in the form: \MSC code \sep code
%% or \MSC[2008] code \sep code (2000 is the default)
Numerical Analysis 
\sep
Wave Propagation
\sep 
Earthquake
\sep
Liquefaction
\sep
Constitutive modeling
\end{keyword}

\end{frontmatter}

%\linenumbers

%%
%% Start line numbering here if you want
%%
% \linenumbers

%% main text


%===============================================

\section{INTRODUCTION}
\label{}
Performance-Based Design (PBD) of geotechnical structures is gaining popularity in professional practice and is fostering research in the academic community.
%
PBD relies heavily on modeling and simulation tools.
%
One of the challenges in PBD of geotechnical and geophysical problems is analysis of dynamic transient phenomena in fluid-saturated porous media and in particular modeling and simulations of seismic wave propagation.
%
This subject, that can be related to liquefaction, continues to challenge engineering research and practice.
%
In addition, lateral movement of sloping ground as a result of pore pressure generation poses other challenges in PBD.


Pore pressure generation phenomenon commonly occurs in loose to medium dense sands which are fully saturated and may induce two related phenomena: the flow liquefaction which leads to flow slides, and the cyclic mobility which leads to lateral spreads.
%
Flow liquefaction occurs when shear stresses required for static equilibrium exceed residual shear strength , as opposed to liquefaction which implies a zero effective stress.
%
In this phenomena earthquake brings soil to point of instability and deformations are driven by static stresses.
%
The associated failure occurs rapidly with little warning and produces large deformations.
%
On the other hand cyclic mobility occurs when shear stresses required for static equilibrium are less than residual shear strength.
%
These deformations are driven by dynamic stresses, occur incrementally and could be large or small \cite{Kramer96}.
%
Accumulation of permanent deformations, degradation of soil moduli, increase of hysteretic damping and change of soil fabric as a function of imposed cyclic shear strains require advanced models and implementations which are tricky but doable.


In order to systematically investigate this subject it is necessary to predict the generation, redistribution, and dissipation of excess pore pressures during and after earthquake shaking and their impact on the transmitted waves. 
%
A fundamental approach requires a dynamic coupled stress-field analysis. 
%
Fully coupled transient response of solid -- pore fluid interaction and constitutive behavior of the soil skeleton play equally important roles in successful numerical simulation of response in saturated granular soil medium.
%
The mechanical model of this interaction when combined with suitable constitutive description of the solid phase and with efficient, discrete, computation procedures, allows most transient and static problems involving deformations to be properly modeled and accurately simulated.


An advanced numerical simulation tool has been developed and utilized for this purpose.
%
The issues of interests here are the capabilities of the advanced constitutive model and the rigorous framework that is used for modeling the fully coupled solid skeleton-pore fluid interaction.
%
Many constitutive models with different levels of complexity have been formulated to describe the response of soils in cyclic loading, e.g. \cite{Martin75,Finn77,Zienkiewicz85,Prevost85_model,Iai92,Manzari97}. 
%
The bounding surface models have generally proved  efficient and successful in simulations of cyclic loading of sands \cite{Dafalias94}. 
%
Recently, greater consideration has been given to the role of the progressive decay in soil stiffness with increasing pore pressure, accumulation of deformation, stress-dilatancy and hysteretic loops in formulating constitutive models for liquefiable soils \cite{Elgamal03,Dafalias04DM}, among others. 


The bounding surface constitutive model in this study predicts with accuracy the soil response of non-cohesive soils during loading under different paths of monotonic and cyclic loading, drained and undrained conditions, and for various soil densities, stress levels and loading conditions. 
%
From a practical point of view, it is most important probably that for all aforementioned conditions, a single soil-specific set of constants is needed.
%
Besides the superiority of this constitutive model in modeling of displacement proportional damping, the other attractive feature of the present framework is proper modeling of velocity proportional damping. This is done by taking into account the interaction of pore fluid and solid skeleton. 
%
In this way the present formulation models in a realistic way the physical damping, which dissipates energy during any wave propagation.


This paper summarizes the key modeling features of our modeling and simulation tool, including the soil stress-strain behavior and the finite element formulation.
%
A detailed calibration of the employed constitutive model for three different sands has been presented which can be used as a guideline for future applications of the model. 
%
Using these model parameters, performance of the constitutive model in reproducing the material response for different types and conditions of loading, as well as different densities and stress states is illustrated.
%
This is followed by verification of the implemented fully coupled element in modeling of a complex step loading on a saturated elastic medium. 



The rest of the paper deals with modeling and simulation of earthquake-induced soil liquefaction and its effects on performance of layered soil systems.
%
Specifically the seismic response of a uniform soil column with relatively dense sand is compared with that of a layered soil column with a liquefiable loose layer in depth.
%
This liquefiable layer could operate as an isolating layer, causing reduction of accelerations in shallow layers. 
%
In addition an example for illustration of the possible effects of a loose interlayer as well as the magnitude of the base acceleration in a relatively dense sloping soil column on the resulting lateral displacements and details of the wave propagation have been numerically studied and discussed. 
%
The presented simulations rely on verified formulation and implementation of behavior of fully coupled porous media and on validated constitutive material modeling. 
%
The numerical simulations presented in this paper illustrate the ability of advanced computational geomechanics tools in providing valuable detailed informations.


%===============================================


\section{MATHEMATICAL FORMULATION}
\label{Formulation}


The material behavior of soil skeleton is interrelated to the pore fluid
pressures. The behavior of pore fluid is assumed to be elastic and thus all
the material nonlinearity is concentrated in the soil skeleton. The soil
behavior (mix of soil skeleton and pore fluid) can thus be described using
single--phase constitutive analysis approach for skeleton combined with the
full coupling with pore fluid. These are the two major parts of the formulation 
in the present study. 


\subsection{Constitutive Model}
\label{}
SANISAND constitutive model is used here for modeling of soil response. SANISAND is the name used for a family of {\small S}imple {\small ANI}sotropic {\small SAND} constitutive models within the frameworks of critical state soil mechanics and bounding surface plasticity. 
%
Manzari and Dafalias \cite{Manzari97} constructed a simple stress-ratio controlled constitutive model for sand in a logical sequence of simple steps. The model is fully compatible with critical state soil mechanics principles by rendering the slope of the dilatancy stress ratio  (also known as the phase transformation line), a function of the state parameter $\psi$, such that at critical state where $\psi=0$, the dilatancy stress ratio coincides with the critical state failure stress ratio. In addition, softening of dense samples is modeled within a collapsing peak-stress ratio bounding surface formulation. The peak stress ratio is again made function of $\psi$ such that at critical state where $\psi=0$ it becomes the critical state stress ratio following an original suggestion by Wood et al. \cite{Wood94}.

 
The bounding surface feature enables simulation of reverse and cyclic loading response simulation. 
%
Dafalias and Manzari \cite{Dafalias04DM} extended the model to account for the fabric changes during the dilatant phase of deformation on the subsequent contractant response upon load increment reversals. 
%
An additional mechanism was also proposed by Dafalias et al. \cite{Dafalias04DPL} for the inherent anisotropy effects. 
%
Using the concept of Limiting Compression Curve \cite{Pestana95} and a proper closed yield surface, Taiebat and Dafalias \cite{Taiebat08SANISAND} eliminated the inability of the previous versions of the model to induce plastic deformation under constant stress-ratio loading, which is especially important at high confining pressures causing grain crushing.


In the present paper the focus is on wave propagation in granular media. In order to involve smaller number of model parameters and for simplicity the version of the SANISAND model with fabric change effects \cite{Dafalias04DM} has been considered as the constitutive model for the soil. The inherent anisotropy \cite{Dafalias04DPL} and the plastic strains under constant-stress ratios \cite{Taiebat08SANISAND} have not be accounted for in the present work. 
%
An outline of the constitutive model in its generalized form for multiaxial stress space is given in the Appendix \ref{Appendix_ConstitutiveModel} for quick reference in the following sections.


\subsection{Fully Coupled Solid Skeleton-Pore Fluid Interaction}
\label{}


The mechanical model of the interaction between solid skeleton and the pore
fluid when combined with suitable constitutive description of the solid
phase and with efficient, discrete, computation procedures, allows one to solve most
transient and static problems involving deformations. The
modeling framework described here is appropriate for saturated porous media,
based on the concepts originally outlined by Biot \cite{Biot56}. 
%
For modeling of the fully coupled solid--pore fluid interaction, Zienkiewicz and Shiomi \cite{Zienkiewicz84} proposed three general continuum approaches, namely (a) $u-p$, (b) $u-U$, and (c) $u-p-U$ formulations (u: solid displacements, p: pore pressure, and U: fluid displacements). 


In this study the coupled dynamic field equations with $u-p-U$ formulation have been used to determine pore fluid and soil skeleton responses. 
%
This formulation uses pore fluid pressure as an additional (dependent) unknown field to stabilize the solution of the coupled system. 
%
The pore fluid pressures have been connected to (dependent on) the displacements of pore fluid so that with known volumetric compressibility of the pore fluid, pressure can be calculated. 
%
Despite it's power, this formulation has rarely been implemented into finite element codes. 
%
The formulation takes into account velocity proportional damping (usually called viscous damping) by proper modeling of coupling of pore fluid and solid skeleton, while the displacement proportional damping is appropriately modeled using elasto--plasticity with the powerful material model chosen. 
%
No additional (and artificial) Rayleigh damping has been used in the FEM model.
%
The main components of the $u-p-U$ formulation for porous media 
are outlined in Appendix  \ref{Appendix_FiniteElementFormulation} 
for quick reference.
%
Detailed description of the $u-p-U$ formulation, finite element discretization and
time integration are presented in \cite{Jeremic08IJNAMG}. 


%===============================================


\section{NUMERICAL TOOL}
\label{}


\subsection{Program Implementation}
\label{}


In order to study the dynamic response of saturated soil systems as a boundary-value problem, an open-source three dimensional numerical simulation tool has been developed. 
%
Parts of OpenSees framework \cite{McKenna97} have been used to connect the finite element domain. 
%
In particular, Finite Element Model classes from OpenSees (namely, class abstractions for Node, Element, Constraint, Load and Domain) have been used to describe the finite element model and to store the results of the analysis performed on the model. 
%
In addition, Analysis classes were used to drive the global level finite element analysis, i.e., to form and solve the global system of equations. 
%
As for the Geomechanics modules, a number of elements, algorithms and material models from UC Davis Computational Geomechanics toolset have been used. 
%
In particular, set of NewTemplate3Dep \cite{Cheng06} numerical libraries have been used for constitutive level integrations, nDarray numerical libraries \cite{Jeremic98nDarray} were used to handle vector, matrix and tensor manipulations, while FEMtools libraries  ($u-p-U$ element) were used in implementation of the coupled solid-fluid interaction at the finite element level. 
%
Finally, solvers from the uMfPACK set of libraries \cite{Davis97} were used to solve the nonsymmetric global (finite element level) system of equations.


The SANISAND constitutive model and the fully coupled $u-p-U$ element are implemented in the above mentioned numerical simulation tool \cite{Cheng06,Taiebat08thesis,Jeremic08IJNAMG}. 
%
The constitutive model is integrated using a refined explicit integration method with automatic error control for yield surface drift based on the work by Sloan et al. \cite{sloan01}.
%
In order to develop integration of dynamic finite element equation in the time domain, the weak form of the finite element formulation is rewritten in a residual form \cite{Argyris91} and the resulting set of residual (nonlinear) dynamic equations is solved using the Hilber--Hughes--Taylor (HHT) $\alpha$--method \cite{Hilber77,Hughes78a,Hughes78b}.
%
The developed finite element model for this formulation uses eight node brick elements. 
%
Because the pore fluid is compressible, there are no problems with locking, particularly if good equation solvers are used, therefor there is no need for lower order of interpolation for pore fluid pressures.
%
We use $u-p-U$ formulation where additional unknown for fluid is used (pore pressure in addition to fluid displacements) and that helps in stabilizing the system. 
%
Pore fluid is compressible but it is many times (orders of magnitude) more stiff than soil skeleton.


All of the above libraries and implementations are available either through their original developers or through second author's website (http://geomechanics.ucdavis.edu).


\subsection{Verification and Validation}
\label{}


Prediction of mechanical behavior comprises use of computational model to foretell the state of a physical system under conditions for which the computational model has not been validated \cite{Oberkampf02}. Confidence in predictions relies heavily on proper Verification and Validation process.


{\bf Verification} is the process of determining that a model implementation accurately represents the developer's conceptual description and specification. It is a Mathematics issue. Verification provides evidence 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.


{\bf Validation} is the process of determining the degree to which a model is accurate representation of the real world from the perspective of the intended uses of the model. It is a Physics issue. {Validation provides evidence that an appropriate model is used.} Validation serves two goals, (a) tactical goal in identification and minimization of uncertainties and errors in the computational model, and (b) strategic goal in increasing confidence in the quantitative predictive capability of the computational model.


Verification and Validation procedures are hence the primary means of assessing accuracy in modeling and computational simulations. 
%
Fig. \ref{Fig-VV} depicts relationships between verification, validation and the computational solution.



In order to verify the $u-p-U$ formulation, a number of closed form or very accurate solutions should be used. To this end and also for illustration of the performance of the formulation and versatility of the numerical implementation, several elastic one-dimensional problems have been studied including
%
drilling of a well, 
the case of a spherical cavity, 
consolidation of a soil layer, 
line injection of fluid in a reservoir, and 
shock wave propagation in saturated porous medium.
%
The verification process for shock wave propagation in porous medium which is the most important and difficult example among these five has been presented here. 
%
Validation has been carried out on an extensive set of physical tests on Toyoura, Nevada, and Sacramento river sands and has been presented along with detailed explanation on calibration of the constitutive model parameters.



\subsubsection{Verification: $u-p-U$ Element}
\label{ShockWave}


The analytic solution developed by Gajo \cite{Gajo95b} and Gajo and Mongiovi \cite{Gajo95} for 1D shock wave propagation in elastic porous medium was used in the verification study. A model was developed consisting of 1000 eight node brick elements, with boundary conditions that mimic 1D behavior (Fig. \ref{ShockWaveMesh}). In particular, no displacement of solid ($u_x = 0 , u_y = 0$) and fluid ($U_x = 0 , U_y = 0$) in $x$ and $y$ directions, i.e. in horizontal directions, was allowed along the height of the model. Bottom nodes have full fixity for solid ($u_i = 0$) and fluid ($U_i=0$) displacements while all the nodes above base are free to move in $z$ direction for both solid and fluid. Pore fluid pressures are free to develop along the model. Loads to the model consist of a unit step function (Heaviside) applied as (compressive) displacements to both solid and fluid phases of the model, with an amplitude of $0.001$~cm. The $u-p-U$ model dynamic system of equations was integrated using the HHT algorithm.


Table~\ref{Table_shock_wave} presents relevant parameters for this verification. Two set of permeability of material were used in the verification. The first model had permeability $k=10^{-6}$~cm/s which creates very high coupling between porous solid and pore fluid. The second model had permeability $k=10^{-2}$~cm/s which creates a low coupling between porous solid and pore fluid. Good agreement was obtained between the numerical simulations and the analytical solution as presented in Figure~\ref{DynamicVerification}.


\subsubsection{Validation: SANISAND Model}


Material model validation can be performed by comparing experimental (physical) results and numerical (constitutive) simulations for different type of sands in different densities, confining pressures, loading paths, and drainage conditions. 
%
The SANISAND constitutive model has been specifically tested for its ability to reproduce a series of monotonic and cyclic tests on Toyoura sand \cite{Verdugo96,Pradhan89_test}, Nevada sand \cite{Arulmoli92}, and Sacramento river sand \cite{Lee67} in a wide range of relative densities and confining pressures and in different drainage conditions. 
%
The 14 material parameters required by the model for these sands are listed in Table \ref{Table-SANISAND} divided in different groups according to the particular role they play. 
%
The material parameters can be selected mainly from standard types of laboratory tests. 
%
Here the calibration has been carried out using drained and undrained triaxial compression
and extension tests at different values of initial void ratio and confining pressure in order to calibrate different features of hardening/softening and dilatancy/contractancy of the model.
%
A step-by-step calibration process for the model constants for Toyoura sand is shown in Figs. \ref{Fig-Toyoura-Calibration-G0}--\ref{Fig-Toyoura-Calibration-Fabric}. 


The $G_{0}$ parameter which defines the elastic shear modulus of sand in
Eq. (\ref{Eq-GK})$_1$ can be calibrated using the stress-strain
curves or from the elastic wave propagation tests in the field or
laboratory. In this paper, a good estimation of $G_{0}$ for monotonic
shearing in Toyoura sand is obtained by fitting the triaxial versions of 
Eqs. (\ref {Eq-Elasticity})$_1$ and (\ref{Eq-GK})$_1$ to the initial stages 
of the deviatoric stress-strain ($\varepsilon_{q}-q$) data of triaxial drained tests in Fig. \ref{Fig-Toyoura-Calibration-G0}. 


The critical state parameters consist of $M_{c,e}$, the critical
stress ratio in triaxial compression and extension, and the parameters $\left(
e_{0}\right)$, $\lambda$, and $\xi$ of Equation $e_{c}=e_{0}-\lambda \left(
p_{c}/p_{at}\right) ^{\xi }$ that defines the position of the CSL in the void
ratio-mean effective stress \cite{Li98}. Calibration of these parameters requires
monotonic tests that approach critical state. Results of the drained and undrained triaxial
compression tests on Toyoura sand are presented in Fig. \ref{Fig-Toyoura-Calibration-CSL} for fitting these parameters.


These bounding and dilatancy surfaces in the SANISAND model are generalization of the 
bounding (or peak), dilatancy (or phase transformation) lines in triaxial space, expressed analytically by 
$\eta=M\exp(-n^b\psi)$ and $\eta=M\exp(n^d\psi)$, respectively. 
%
The parameters $n^{b}$ and $n^{d}$ can be determined by evaluating these equations at the peak and phase transformation states. The calibration process is shown in Fig. \ref{Fig-Toyoura-Calibration-BoundingDilatancy}(a,b) where 
$\psi^{b}$ and $\eta^{b}$ are the values of $\psi$ and $\eta$ at the peak stress-ratio state, 
and 
$\psi^{d}$ and $\eta^{d}$ are the values of $\psi$ and $\eta$ at the phase transformation state. 
%
Note that the phase transformation state corresponds
to the peak of the volumetric strain $\varepsilon_{v}$ in drained tests and
the peak of the excess pore pressure $u-u_{0}$ in undrained tests. 
%
Direct estimation of the dilatancy parameter $A_{0}$ requires good
quality stress-dilatancy data, where $A_{0}$ can then be calibrated based 
on $\varepsilon_{v}-\varepsilon_{q}$ curves. Here the parameters 
$A_0$ (related to dilatancy), and $h_{0}$ and $c_{h}$ (related to the effect of 
distance from bounding surface) are estimated by trial-and-error as shown in Fig. \ref{Fig-Toyoura-Calibration-BoundingDilatancy}(c-f).


Finally fabric-dilatancy constants $z_{max}$ and $c_z$ require trial-and-error 
fitting of loading--unloading reverse loading, or cyclic 
data, preferably undrained, where $\eta$ must exceed $M^d$ in the 
process so that evolution of $z$ is activated.
%
Fig. \ref{Fig-Toyoura-Calibration-Fabric} shows calibration of the $z_{max}$ parameter using the loading--unloading undrained triaxial tests on Toyoura sand.


Figs. \ref{Fig-Toyoura-CIDC} and \ref{Fig-Toyoura-CIUC} compare the experimental data \cite{Verdugo96} and model simulations for drained and undrained triaxial compression tests followed by unloading on isotropically consolidated samples of Toyoura sand. These experiments/simulations cover a wide range of confining pressures and densities, and show variations of response from highly dilatant to highly contractant. This variety of response depends on the combination of density and confining pressure that is reflected in the value of state parameter $\psi$. 


Figs. \ref{Fig-Toyoura-CYDC-loose} and \ref{Fig-Toyoura-CYDC-dense} compare the experimental data \cite{Pradhan89_test} and model simulations for drained constant-p cyclic triaxial tests on isotropically consolidated
loose and dense samples of Toyoura sand. In these tests the amplitude of shear strain ($\gamma=\epsilon_{a}-\epsilon_{r}$ as defined in \cite{Pradhan89_test}) has been increased with cyclic loading. Accumulation of compressive volumetric strain with cyclic loading can be observed in Fig. \ref{Fig-Toyoura-CYDC-loose}, and clearly when the stress ratio exceeds a certain value, which varies by the state, the specimens begins to dilate. As it is expected the dilatancy feature in more pronounced for the dense sample in Fig. \ref{Fig-Toyoura-CYDC-dense}.


Experimental results \cite{Arulmoli92} and numerical simulation for constant-p triaxial compression tests on Nevada sand in different densities and confining pressures are presented in Fig. \ref{Fig-Nevada-CIDC}. In addition, the stress path and shear response of an undrained triaxial test on anisotropically consolidated Nevada sand, at relative density of about $40\%$, are presented in Fig. \ref{Fig-Nevada-CYUC}.


The last set of simulations for the purpose of model validation are presented in Fig. \ref{Fig-Sacramento-CIDC} where experimental results of drained triaxial compression tests on Sacramento river sand \cite{Lee67} have been compared with the model simulations. The experiments/simulations are on both loose and dense samples of Sacramento river sand in different initial confining pressures. It should be noted that the fabric dilatancy parameters, needed in correct simulation loading-unloading cycles, have not been calibrated for Sacramento river sand as the authors could not find enough related experimental data in the literature.

 
Close match of physical testing data with constitutive predictions represents a satisfactory validation of the material models. This validation with previous verification gives confidence in simulations of real, prototype behavior.


%===============================================


\section{NUMERICAL SIMULATION EXAMPLES}
\label{NumericalSimulationExamples}
Liquefaction of level and sloping ground represents a common situation during earthquakes. It is therefore of interest to use the developed numerical tool for simulating the earthquake wave propagation is such cases.


This section presents results of numerical simulation for different cases of a 1D site response, using the available 3D elements in the developed numerical simulation tool. Two general examples have been studied. The first example presents results of a study on  the isolating effect of a liquefied sand layer in propagation of seismic waves. In the second example the seismic induced shear deformation of a gentle slope in presence or absence of liquefiable layer has been studied.


It should be noted that the simulations are based on the small deformation assumption, which might introduce large error in strain tensor (by neglecting quadratic portion of displacement derivatives)  when soil undergoes large deformations.


\subsection{Example 1: Isolation of Ground Motions in Level Ground}
\label{Example1}


A 10m vertical column of saturated sand consisting of twenty $u-p-U$ eight-node brick elements was considered subjected to an earthquake shaking at the base. Two particular cases have been studied in this example.
%
In the first case, all 20 elements of the soil column are assumed to be at medium dense state at the beginning of analysis ($e_{in}=0.80$). In the second case, the deepest two meters of the soil column are assumed to be at looser initial states compared to the first case. In particular two of the deepest elements are assumed to be at $e_{in}=0.95$ (loose) and the next two elements at $e_{in}=0.875$ (medium loose), just to make a smoother transition between the loose and dense elements. The rest of the soil column, i.e. the upper 16 elements (8 meters) are at the same density as the first case ($e=0.80$). A schematic illustration of these two cases is presented in Fig. \ref{Fig-Isolation}.


The elements are labeled from E01 (bottom) to E20 (surface) in Fig. \ref{Fig-Isolation}. The boundary conditions are such that the vertical displacement degrees of freedom (DOF) of the soil and water at the base ($z=0$) are fixed, while the pore pressure DOFs are free. The soil and water vertical displacement DOFs at the surface ($z=10m$) are free to simulate the upward drainage. The pore pressure DOFs are fixed at the surface. On the sides, soil skeleton and water are prevented from moving in the $y$ direction while movements in $x$ and $z$ directions are free. 
%
It is emphasized that the displacements of soil skeleton and pore fluid are different. 
%
In order to simulate the 1D behavior, all of the related DOFs of the nodes at each depth are connected in a master--slave fashion, i.e. the nodes at the same level have the same $u_i$ (i=x,y,z). 
%
Same thing is true for $p$ and $U_i$. 
%
The values of $u_i$ and $U_i$ however may be different at a node, allowing for relative movement of  soil and fluid.
%
%This is for each of the phases, that is $u_i$ are same and $U_i$ are same among themselves, but soil and fluid move differently.
%
The soil is assumed to be Toyoura sand and is characterized by the validated SANISAND model explained in previous sections with the model parameters given in Table \ref{Table-SANISAND}.  The other parameters related to the boundary value problem are given in Table \ref{Table-BVP}. The permeability is assumed to be isotropic and equal to $k=5.0{\times}10^{-4}$~m/s.


A self-weight analysis was performed before the base excitation. The resulting fluid hydrostatic pressures and soil stress states along the soil column serve as initial conditions for the subsequent dynamic analysis. An input acceleration time history (Fig. \ref{Fig-Isolation}), taken from the recorded horizontal acceleration of Model No.1 of VELACS project at RPI \cite{Arulanandan93}, was considered in which $a_{max}=0.2g$ and main shaking lasts for about 12 seconds.


In addition to the [physical] velocity and displacement proportional damping from the element and the material model, respectively, a small amount of numerical damping has been introduced through constants of the HHT algorithm in order to stabilize the dynamic time stepping. 
%
No Rayleigh damping has been used.


Results of the analysis during and after the \emph{shaking phase} for both cases (i.e. the uniform and the layered soil profiles) are presented in Figs. \ref{Fig-Isolation-stresspath}--\ref{Fig-Isolation-accel}. 
%
In all these figures the plots on the left and right correspond to the cases of \emph{uniform} and \emph{layered} soil profiles, respectively. 
%
In particular, Figs. \ref{Fig-Isolation-stresspath} and \ref{Fig-Isolation-stressstrain} show variations of the shear stress $\sigma_{xz}$ vs vertical effective stress $\sigma_{z}$ and shear strain $\gamma$, respectively.
%
Figs. \ref{Fig-Isolation-voidratio} and \ref{Fig-Isolation-accel} show time histories of void ratio $e$ and horizontal ($x$) component of acceleration in solid part of the mixture $a_{u,x}$, respectively.
%
Note that different rows of plots in Figs. \ref{Fig-Isolation-stresspath}, \ref{Fig-Isolation-stressstrain}, and \ref{Fig-Isolation-voidratio} correspond to different elements. Similarly, different rows in Fig. \ref{Fig-Isolation-accel} correspond to different elevations of the nodes.


Fig. \ref{Fig-Isolation-stresspath}a shows the typical mechanism of cyclic decrease in vertical effective stress for the case of the uniform soil column. Signs of the so--called butterfly loops in the effective stress path can be observed partially from early stages of shaking at the upper layer as these layers are in lower confining pressures compared to the deeper layers and thus are more dense than critical in the critical state soil mechanics terminology, i.e. with less contractive tendency. In later stages of shaking, when the confining pressures are reduced to smaller values, the butterfly shapes of the stress paths are more pronounced and almost all depths of the soil column experience the cyclic mobility mechanism of liquefaction. In this stage the soil layers momentarily experience small values of effective stress followed by a dilative response in which the elements recover their strength in cycles of loading due to the their denser than critical state. As a result, the seismic-induced shear waves propagate all the way to the surface of the soil column. The response can be observed as an average degradation of stiffness and accumulation of shear strains in all levels of the soil column as shown in Fig. \ref{Fig-Isolation-stressstrain}a. 


In case of the layered soil profile, however, Fig. \ref{Fig-Isolation-stresspath}b shows that the first cycle of shaking degrades the vertical effective stress in the loose layer of element 1, E01, to very small values in a flow liquefaction mechanism. The shear wave in this first cycle of shaking propagates to the surface and causes some reduction in vertical effective stress in upper soil layers. Because element 1, E01, is in loose states and has a strong contractive tendency it does not recover its strength in cyclic loading and therefore it remains in this liquefied state after the first cycle with negligible shear resistance. Therefore this layer acts as an isolating layer and no shear stress can be transmitted to the upper layers after the first cycle. This mechanism can also be observed in the shear stress -- shear strain plot of Fig. \ref{Fig-Isolation-stressstrain}b as the first cycle completely liquefies the E01 and after this cycle no shear stress is propagated to the upper layers. Of course the very loose element E01 shows considerable amount of shear strain due to liquefaction. As it is shown in Fig. \ref{Fig-Isolation-stresspath}b in upper layers after the first cycle of shaking no reduction in vertical effective stress occurs \emph{due to shearing} (induced by shaking) in these layers, however, one can observe that the vertical effective stress keeps decreasing in upper layers during (and even for a while after) shaking. This is because the dissipated pore pressure in deeper layers travels upward and reduces the effective stress in upper layers. The effective stress will then start to recover as the excess pore pressure dissipates from these layers. 


Fig. \ref{Fig-Isolation-voidratio} shows the redistribution of void ratio during and long after the shaking. 
The horizontal dashed line in each plot shows $e_{in,shaking}$ that is the value of void ratio for the corresponding layer at the beginning of shaking (at the end of self-weight analysis).
%
For the uniform soil column Fig. \ref{Fig-Isolation-voidratio}a shows a general trend of consolidation in all of the soil layers. This consolidation starts later in upper layers because of the water pumped from the lower layers, which is to be dissipated by time. 
%
In the layered soil column (Fig. \ref{Fig-Isolation-voidratio}b) the loose layer in element 1, E01, shows considerable consolidation due to shaking. The dissipated volume of water from this layer is pumped to the upper layers and result is some dilation (increase in void ratio) in the these layers for some time after the end of shaking. By the time this pumped volume of water in the upper layers dissipates all the layers show a consolidation response. 
%
In other words in the layered soil column pumping of water causes initial loosening of soil in the upper layers, which is then followed by densification. 


Finally Fig. \ref{Fig-Isolation-accel} shows propagation of the horizontal acceleration through the soil column during shaking. In the uniform soil column Fig. \ref{Fig-Isolation-accel}a shows that the base acceleration is transmitted to the surface of the soil column. It shows some spikes of acceleration with large amplitude at the surface of the soil resulting from the dilative phases of soil response during the shaking phase. In the layered soil column, however, Fig. \ref{Fig-Isolation-accel}b shows that only the first cycle of shaking is transmitted to the surface (although with reduced amplitude)  and subsequently the upper layers are isolated because of the fully liquefied state of the loose layer in element 1, E01. The transmitted accelerations to the surface of the soil column after the first cycle of shaking have very small amplitudes. 


\subsection{Example 2: Earthquake-Induced Shear Deformations in Gentle Slopes}
\label{Example2}


A 10m gently inclined ($1.5^{\circ}$) column of saturated sand consisting of twenty $u-p-U$ eight-node brick elements was considered subjected to a base shaking. In terms of initial densities two particular cases have been studied in this example. In the first, the soil column is assumed to be uniform with $e=0.7$ (relatively dense, $Dr\simeq75\%$), while in the other case a similar soil column is considered which includes a relatively loose layer with $e=0.9$ ($Dr\simeq25\%$) at depth 2.5m to 3.5m (Fig. \ref{Fig-mesh-accel}). 


The boundary condition are exactly same as those presented in Example 1.
%
The soil parameters and the other parameters related to the boundary value problem are also same as those in Example 1, and are given in Tables \ref{Table-SANISAND} (Toyoura sand) and \ref{Table-BVP}. The permeability is assumed to be isotropic and equal to $k=10^{-3}$~m/s.


A static application of gravity analysis was performed before the base excitation. The resulting fluid hydrostatic pressures and soil stress states along the soil column serve as initial conditions for the subsequent dynamic
analysis. The input acceleration time history (Fig. \ref{Fig-mesh-accel}) is considered with $a_{max}=0.2g$ and shaking lasts for about 8 seconds. 
%
The following three cases have been studied in this paper:%
\begin{description}
\item [I.] Uniform soil profile with $a_{max}=0.2g$
\item [II.] Layered soil profile with $a_{max}=0.2g$
\item [III.] Layered soil profile with $a_{max}=0.4g$
\end{description}


The predicted behavior for cases I, II and III in terms of the vertical effective stress, shear stress, shear strain, excess pore pressure ratio ($R_u$) and lateral displacement are presented and discussed in the rest of this section. 


The response is shown for selected elements along the soil profile (for positions of the elements refer to Fig. \ref{Fig-mesh-accel}).
%
In particular Figs. \ref{Fig-case-I}a, \ref{Fig-case-II}a and \ref{Fig-case-III}a show changes in the shear stress and vertical effective stress. Figs. \ref{Fig-case-I}b, \ref{Fig-case-II}b and \ref{Fig-case-III}b show changes in the shear stress and shear strain. Finally Figs. \ref{Fig-case-I}c, \ref{Fig-case-II}c and \ref{Fig-case-III}c show times histories of the excess pore pressure ratio, $R_u$.  Contours of excess pore pressure ($R_u$), shear strain, and lateral displacements in all depths vs. time are presented in parts d, e and f of Fig. \ref{Fig-case-I}, \ref{Fig-case-II} and \ref{Fig-case-III} for the above three cases.
%
Note that $R_u=u_e/(\sigma_v)_{in}$ where $u_e$ is the excess pore pressure and $(\sigma_v)_{in}$ is the initial effective vertical stress. Clearly $R_u$ cannot be larger than $1$ since $u_e\leq(\sigma_v)_{in}$ (in cohesionless soils the stress state cannot go to the tensile side, in which $u_e>(\sigma_v)_{in}$).


Typical mechanism of response observed in different elements during the shaking is explained in the following. 
%
Generally at the early stages of loading in a cycle the element shows a contractive response in shearing which is associated with an increase in excess pore pressure and a decrease in vertical effective stress and the stiffness (in particular the shear stiffness). As loading continues to some larger stress-ratios the mechanism of response in the element might change to a dilative regime. This means that as shearing continues, the excess pore pressure decreases and the vertical effective stress and the soil stiffness increase. In other words soil regains its stiffness during the dilative phase. Upon and after reversal of shear stress increment again the sample shows a contractive regime with an increase in excess pore pressure and a decrease in vertical effective stress and soil stiffness. As the absolute value of  stress-ratio exceeds a certain level (phase transformation) in the reverse loading (in negative stress-ratios) the element again shows the dilative response. The final stress-ratio increment reversal in the loading cycle is again associated with a contractive response. 
%
The above mechanisms typically happens in all cycles. The overall contractive response could be stronger or weaker than the dilative response depending on the state of soil, i.e. its void ratio and effective stress. 
%
For looser samples or at higher confining pressures the contractive response is more pronounced and therefore shearing at such states shows increase in excess pore pressure and decrease in vertical effective stress and stiffness. In denser samples or at lower confining pressures the response is the opposite. 


Fig. \ref{Fig-case-I} shows detailed response of the soil profile in case I. 
%
Based on the previous general explanation the response starts with a contractive regime which results in decreasing the vertical effective stress and the stiffness, and increasing the excess pore pressure ratio. 
%
Since a relatively dense state is assumed for this analysis ($e=0.7$, $D_r\simeq75\%$) the soil has a general tendency of dilation as the  stress-ratio exceeds the dilation line (phase transformation). On the other hand, the inclination of the soil column poses asymmetric horizontal shear stresses (toward the direction of inclination) during cycles of shaking. 
%
This essentially makes the sample more prone to dilation in the parts of shaking toward the the inclination direction (positive directions of shear stress and strain in this case) and the soil regains its stiffness and strength ($\sigma_v$) in the dilative parts of the loading cycle. 
%
The dilative tendency of response besides the offset shear stress in cyclic loading lead to asymmetric butterfly loops in stress space during the dilative and the subsequent intense contractive phases of response. 
%
The asymmetric horizontal shear stresses also results in asymmetric shear strains which accumulate toward the inclination of the soil profile and create permanent horizontal displacements in this direction. 
%
The dilative phases lead to instantaneous drops in excess pore pressure ratio $R_{u}$. These drops are more pronounced in upper elements which have lower confining pressures and therefore have stronger dilative tendency.


Fig. \ref{Fig-case-II} shows response of the soil profile in case II. In this case presence of the relatively loose layer at depth 2.5-3.5m ($e=0.9$, $D_r\simeq25\%$) affects the response of the system. This layer has stronger contractive response in cycles of loading. Such tendency causes faster increase in excess pore pressure and decrease in vertical effective stress and stiffness. Degradation of stiffness leads to faster accumulation of the permanent shear strain in the symmetric loops of stress -- strain. The time histories of surface lateral displacement for these two cases are shows in Fig. \ref{Fig-ux-compare} (presented results of case III in this figure will be discussed later). It can be observed that the permanent surface lateral displacement in case II ($\simeq0.40m$) is almost twice as much as case I  ($\simeq0.19m$) due to presence of the relatively loose layer. The excess pore pressure ratio in the loose layers shows instantaneous spikes to $R_u=1$ and drops mainly to two values of $\simeq0.35$ or $\simeq0.65$ which are related to the end of dilative phases in positive and negative shear stresses, respectively. The expelled pore water from this layer during shaking is pumped to the upper layers and shows a considerable increase in the excess pore pressure ratio of element 19 (E19) after the end of shaking. 


Fig. \ref{Fig-case-III} shows response of the soil profile in case III. The geometry and parameters in case III are the same as those in case II. The only difference is the intensity of the input base acceleration which is 0.4g in case III. The interesting observation in the simulation results is that although the magnitude of the base acceleration in case III is twice as much as case II, the calculated residual shear strains are not considerably different in the two cases. 
%
This could be attributed to the fact that in case III the larger amplitude of acceleration (and shear stress) soon brings the soil state to the dilative phase during which the vertical effective stress increases and soil regains its strength. Comparing the stress path of element 03 (E03) in cases II and III makes this phenomena more clear. Due to increase of the vertical effective stress during the intense dilation in case III the effective stress state remains larger in this case compared to case II and therefore less degradation is observed in the subsequent cycles of loading.  The strongly pronounced dilative response in all layers of soil profile in case III clearly affects the resulting excess pore pressure ratios (Fig. \ref{Fig-case-III}c) comparing to what was observed in case II (Fig. \ref{Fig-case-II}c). More specifically, case III experiences more intense drops in plots of $R_u$ even to negative values, specially at upper elements in which are even more dilative due to lower surcharge (confining pressure). 


The time histories of the surface lateral displacement, and variations with depth of the permanent lateral displacement for different cases are presented and compared in Fig. \ref{Fig-ux-compare}. Although case III shows larger amplitudes of cyclic motions, its permanent surface lateral displacement is about 20\% smaller than case II. %
That is quite interesting, and make for the importance of interaction of dynamic characteristics of the earthquake and the response of the soil (and its components) as described by the SANISAND constitutive model.
%
Fig. \ref{Fig-ux-compare}b provides more detailed information about the lateral displacements of soil layers in different cases. Due to the aforementioned stronger dilative response in case III, the dense layers adjacent to the loose layer appear to get stronger in case III and show lesser lateral displacements. This could also be partially explained by smaller surface lateral displacement in case III compared to case II. Note that the permanent displacement is equivalent to integration of shear strains in the whole soil column. 


%===============================================


\section{SUMMARY AND CONCLUSIONS}
\label{}


An efficient finite element formulation and a numerical tool for analysis of wave propagation phenomena in fluid--saturated porous media has been developed.
%
The saturated porous medium is modeled as a two--phase, fully coupled system consisting of a porous solid and a fluid phase.
%
The numerical tool includes the $u-p-U$ formulation for fully coupled behavior of soils and the SANISAND model as an advanced elastic--plastic constitutive model for modeling of
stress-strain response of the solid phase.
%
The formulation takes into account velocity proportional damping (usually called viscous damping) by proper modeling of coupling of pore fluid and solid skeleton, while the displacement proportional damping is appropriately modeled by the dissipation mechanism of the constitutive model of elasto--plasticity used.


Numerical results which demonstrate the accuracy and versatility of the formulations and capabilities of the model are presented.
%
The verified and validated models and computational simulations tool are used to predict behavior of
layered soil system during seismic loading.
%
Two general problems are studied.
%
Firstly, a numerical study is conducted on isolating effects of a liquefied sand layer in propagation of seismic waves.
%
Secondly, seismically induced shear deformation of a gentle slope in the presence or absence of liquefiable layer are studied as well.
%
Attention is given to propagation of shear waves through soil layers, as well as to variation of stresses and strains during such propagation.
%
Spatial evolution of liquefaction in loose and medium dense soil layers is also investigated.
%
In the first example it is shown that liquefaction of loose layers in depth prevents transmission of earthquake-induced motions and shear stresses to the upper layers.
%
While prevention of transmission of motions may be considered as a positive feature, the water flux coming from the liquefied bottom layers might force the upper layers towards lower effective stress and possibly instability even in absence of earthquake-induced shear stresses.
%
More importantly, attention needs to be paid to the fact that liquefaction of the deep loose layers could result in large strains in these layers which might result in (rigid block) translation of the upper layers.
%
The second example is used to present  complexities arising from variation of soil density in a layered system as well as the input base acceleration.
%
For example, it is shown that larger shaking produces smaller lateral spread, because of close interaction of earthquake loading and soil constitutive response (with its solid and fluid components).


The capabilities of the $u-p-U$ element, the SANISAND constitutive model, and the nonlinear dynamic finite element numerical tool are demonstrated in handle complex phenomena. These studies provide new insight into the mechanisms of wave propagation and seismic behavior of saturated soil systems.
%
Fully coupled nonlinear dynamic numerical simulations provide accurate means for detailed investigation on complex loading conditions and geometries and different material
states.
%
Such high fidelity simulations provide detailed input to PBD methods that most existing empirical approaches used in liquefaction modeling fail to provide. 
%
In particular one must emphasize the importance of using a robust and realistic constitutive model such as SANISAND which can account for the various important features of the soil response in relation to its density and loading conditions. 
%
For example it is important that the constitutive model with a single set of parameters can show the dilative or contractive response depending on the level of density and confining pressure in the soil. 
%
This is accomplished by using the state parameter $\psi$.
%
This parameters in the model controls the peak stress ratio, and by consequence the plastic modulus, as well as the dilatancy stress ratio (phase transformation stress ratio), and by consequence the dilatancy.
%
Similarly by using a fabric dilatancy tensor the model has the capability to account for intense contractive response upon reverse loading increments after a dilative phase of loading. This feature has an importsnt effect on the predicted displacement and pore water pressure built up during cyclic loading. Had these constitutive properties not been realistically described, the overall description of the phenomenon would have been much poorer to say the least. 

%===============================================

\SingleSpaced

\clearpage

%======================
% table I
\begin{table}[htdp]
\caption{Simulation parameters used for the shock wave propagation verification problem.}
\begin{center}
\begin{tabular}{lccc}
\hline
\textbf{Parameter} & \textbf{Symbol} & \textbf{Value}\\ 
\hline
Poisson ratio & $\nu$ & 0.3 \\
Young's modulus & E & $1.2 \times 10^6~kN/m^2$ \\
Solid particle bulk modulus & $K_s$ & $3.6 \times 10^7~kN/m^2$ \\
Fluid bulk modulus & $K_f$ & $2.17 \times 10^{6}~kN/m^2$ \\
Solid density & $\rho_s$ & $2700~kg/m^3$ \\
Fluid density & $\rho_f$ & $1000~kg/m^3$ \\
Porosity & n & $0.4$ \\
HHT parameter & $\alpha$ & -0.2 \\
\hline
\end{tabular}
\end{center}
\label{Table_shock_wave}
\end{table}%
%======================
% table II
\begin{table}[htdp]
\caption{Material parameters of the SANISAND constitutive model for three types of sands.}
\begin{center}
\begin{tabular}{lcccc}
\hline
\textbf{Parameter} & \textbf{Symbol} & \textbf{Toyoura} & \textbf{Nevada}
& \textbf{Sacramento} \\ \hline
Elasticity & $G_{0}$ & $125$ & $150$ & $200$ \\ 
& $\nu $ & $0.05$ & $0.05$ & $0.2$ \\ 
CSL & $M$ & $1.25$ & $1.14$ & $1.35$ \\ 
& $c$ & $0.712$ & $0.78$ & $0.65$ \\ 
& $e_{0}$ & $0.934$ & $0.83$ & $0.96$ \\ 
& $\lambda $ & $0.019$ & $0.027$ & $0.028$ \\ 
& $\xi $ & $0.7$ & $0.45$ & $0.7$ \\ 
Dilatancy & $n^{d}$ & $2.1$ & $1.05$ & $2.0$ \\ 
& $A_{0}$ & $0.704$ & $0.81$ & $0.8$ \\ 
Kinematic & $n^{b}$ & $1.25$ & $2.56$ & $1.2$ \\ 
Hardening & $h_{0}$ & $7.05$ & $9.7$ & $5.0$ \\ 
& $c_{h}$ & $0.968$ & $1.02$ & $1.03$ \\ 
Fabric dilatancy & $z_{\max }$ & $2.0$ & $5.0$ & $-$ \\ 
& $c_{z}$ & $600$ & $800$ & $-$ \\ 
\hline
\end{tabular}
\end{center}
\label{Table-SANISAND}
\end{table}%
%======================
% table III
\begin{table}[htdp]
\caption{Parameters used in the boundary value problem simulations (in addition to the material model parameters).}
\begin{center}
\begin{tabular}{lccc}
\hline
\textbf{Parameter} & \textbf{Symbol} & \textbf{Value}\\ 
\hline
Solid particle bulk modulus & $K_s$ & $3.6 \times 10^7~kN/m^2$ \\
Fluid bulk modulus & $K_f$ & $2.17 \times 10^{6}~kN/m^2$ \\
Solid density & $\rho_s$ & $2700~kg/m^3$ \\
Fluid density & $\rho_f$ & $1000~kg/m^3$ \\
HHT parameter & $\alpha$ & -0.2 \\
\hline
\end{tabular}
\end{center}
\label{Table-BVP}
\end{table}%
%======================


\clearpage


%======================
\begin{figure}[pb]
\begin{center}
\includegraphics[angle=-90,width=0.48\textwidth]{VV.eps}
\end{center}
\par
\caption{Schematic description of Verification and Validation procedures and the computational solution  (after \cite{Schlesinger79,Oberkampf02}).}
\label{Fig-VV}
\end{figure}
%======================
\begin{figure}[pb]
\begin{center}
\vspace{-1.0cm}
\includegraphics[angle=-90,width=0.6\textwidth]{ShockWaveMesh.eps}
\end{center}
\par
\vspace{-0.5cm}
\caption{Viscous coupling finite-element analysis of infinite half-space subjected at the surface to a step displacement boundary condition.}
\label{ShockWaveMesh}
\end{figure}
%======================
\begin{figure}[pb]
\begin{center}
\includegraphics[width=0.99\textwidth]{NewMark_gamma0.6_K10-2and10-6cm-THESIS.eps}
\end{center}
\par
\caption{Compressional wave in both solid and fluid, comparison with closed form solution.}
\label{DynamicVerification}
\end{figure}
%======================
\begin{figure}[pb]
\begin{center}
{\normalsize 
\mbox{
\subfigure{\includegraphics[width=0.48\textwidth]{Elasticity.eps}} 
} \\
}
\end{center}
\par
\caption{Calibration of the $G_0$ constant using monotonic triaxial tests on Toyoura sand \cite{Verdugo96}.}%
\label{Fig-Toyoura-Calibration-G0}%
\end{figure}
%======================
\begin{figure}[pb]
\begin{center}
{\normalsize 
\mbox{
\subfigure[]{\includegraphics[width=0.48\textwidth]{CSL-1.eps}}  \quad 
\hspace{-1cm}
\subfigure[]{\includegraphics[width=0.48\textwidth]{CSL-2.eps}} 
} \\
}
\end{center}
\par
\caption{Calibration of the CSL constants using monotonic triaxial tests on Toyoura sand \cite{Verdugo96}. }%
\label{Fig-Toyoura-Calibration-CSL}%
\end{figure}
%======================
\begin{figure}[pb]
\begin{center}
{\normalsize 
\mbox{
\subfigure[]{\includegraphics[width=0.48\textwidth]{nd-nb-1.eps}}  \quad 
\hspace{-1cm}
\subfigure[]{\includegraphics[width=0.48\textwidth]{nd-nb-2.eps}} 
} \\
\mbox{
\subfigure[]{\includegraphics[width=0.48\textwidth]{A0-Toyoura-1.eps}}  \quad 
\hspace{-1cm}
\subfigure[]{\includegraphics[width=0.48\textwidth]{A0-CIUC-Toyoura-1.eps}} 
} \\
\mbox{
\subfigure[]{\includegraphics[width=0.48\textwidth]{h-Toyoura-2.eps}}  \quad 
\hspace{-1cm}
\subfigure[]{\includegraphics[width=0.48\textwidth]{h-interpolation-Toyoura.eps}} 
} \\
}
\end{center}
\par
\caption{Calibration of the bounding and dilatancy surface constants using monotonic triaxial tests on Toyoura sand \cite{Verdugo96}.}%
\label{Fig-Toyoura-Calibration-BoundingDilatancy}%
\end{figure}
%======================
\begin{figure}[pb]
\begin{center}
{\normalsize 
\mbox{
\subfigure{\includegraphics[width=0.48\textwidth]{fabric-Toyoura-1.eps}} 
} \\
}
\end{center}
\par
\caption{Calibration of the $z_{max}$ parameter using the  loading--unloading triaxial tests on Toyoura sand \cite{Verdugo96}. }%
\label{Fig-Toyoura-Calibration-Fabric}%
\end{figure}
%======================
\begin{figure}[pb]
\begin{center}
{\normalsize 
\mbox{
\subfigure{Experiment} \quad 
\hspace{4cm}
\subfigure{Simulation}
} \\ 
\setcounter{subfigure}{0}%
\mbox{
\subfigure[]{\includegraphics[width=0.48\textwidth]{CIDC-1_test.eps}} \quad 
\hspace{-1cm}
\subfigure[]{\includegraphics[width=0.48\textwidth]{CIDC-1_model.eps}}
} \\ 
\mbox{
\subfigure[]{\includegraphics[width=0.48\textwidth]{CIDC-2_test.eps}} \quad 
\hspace{-1cm}
\subfigure[]{\includegraphics[width=0.48\textwidth]{CIDC-2_model.eps}}
}}\\
\end{center}
\par
\caption{
Simulations versus experiments in drained triaxial compression
tests on isotropically consolidated samples of Toyoura sand \cite{Verdugo96}.
}%
\label{Fig-Toyoura-CIDC}%
\end{figure}
%======================
\begin{figure}[pb]
\begin{center}
{\normalsize 
\mbox{
\subfigure{Experiment} \quad 
\hspace{4cm}
\subfigure{Simulation}
} \\ 
\setcounter{subfigure}{0}%
\mbox{
\subfigure[]{\includegraphics[width=0.48\textwidth]{CIUC-1_test.eps}} \quad 
\hspace{-1cm}
\subfigure[]{\includegraphics[width=0.48\textwidth]{CIUC-1_model.eps}}
} \\ 
\mbox{
\subfigure[]{\includegraphics[width=0.48\textwidth]{CIUC-2_test.eps}} \quad 
\hspace{-1cm}
\subfigure[]{\includegraphics[width=0.48\textwidth]{CIUC-2_model.eps}}
}\\
}
\end{center}
\par
\caption{ 
Simulations versus experiments in undrained triaxial compression
tests on isotropically consolidated samples of Toyoura sand \cite{Verdugo96}.
}%
\label{Fig-Toyoura-CIUC}%
\end{figure}
%======================
\begin{figure}[pb]
\begin{center}
{\normalsize 
\mbox{
\subfigure{Experiment} \quad 
\hspace{4cm}
\subfigure{Simulation}
} \\ 
\setcounter{subfigure}{0}%
\mbox{
\subfigure[]{\includegraphics[width=0.48\textwidth]{CIDCp-a_test.eps}} \quad 
\hspace{-1cm}
\subfigure[]{\includegraphics[width=0.48\textwidth]{CIDCp-a_model.eps}}
} \\ 
\mbox{
\subfigure[]{\includegraphics[width=0.48\textwidth]{CIDCp-b_test.eps}} \quad 
\hspace{-1cm}
\subfigure[]{\includegraphics[width=0.48\textwidth]{CIDCp-b_model.eps}}
} \\ 
\mbox{
\subfigure[]{\includegraphics[width=0.48\textwidth]{CIDCp-c_test.eps}} \quad 
\hspace{-1cm}
\subfigure[]{\includegraphics[width=0.48\textwidth]{CIDCp-c_model.eps}}
} \\ 
}
\end{center}
\par
\caption{ 
Simulations versus experiments in constant-p cyclic triaxial 
tests on a relatively loose sample of Toyoura sand \cite{Pradhan89_test}.
}%
\label{Fig-Toyoura-CYDC-loose}%
\end{figure}
%======================
\begin{figure}[pb]
\begin{center}
{\normalsize 
\mbox{
\subfigure{Experiment} \quad 
\hspace{4cm}
\subfigure{Simulation}
} \\ 
\setcounter{subfigure}{0}%
\mbox{
\subfigure[]{\includegraphics[width=0.48\textwidth]{CIDCp_dense-a_test.eps}} \quad 
\hspace{-1cm}
\subfigure[]{\includegraphics[width=0.48\textwidth]{CIDCp_dense-a_model.eps}}
} \\ 
\mbox{
\subfigure[]{\includegraphics[width=0.48\textwidth]{CIDCp_dense-b_test.eps}} \quad 
\hspace{-1cm}
\subfigure[]{\includegraphics[width=0.48\textwidth]{CIDCp_dense-b_model.eps}}
} \\ 
\mbox{
\subfigure[]{\includegraphics[width=0.48\textwidth]{CIDCp_dense-c_test.eps}} \quad 
\hspace{-1cm}
\subfigure[]{\includegraphics[width=0.48\textwidth]{CIDCp_dense-c_model.eps}}
} \\ 
}
\end{center}
\par
\caption{ 
Simulations versus experiments in constant-p cyclic triaxial 
tests on a relatively dense sample of Toyoura sand \cite{Pradhan89_test}.
}%
\label{Fig-Toyoura-CYDC-dense}%
\end{figure}
%======================
\begin{figure}[pb]
\begin{center}
{\normalsize 
\mbox{
\subfigure[]{\includegraphics[width=0.48\textwidth]{CIDC-Nevada-1.eps}} \quad 
\hspace{-1cm}
\subfigure[]{\includegraphics[width=0.48\textwidth]{CIDC-Nevada-2.eps}}
} \\ 
}
\end{center}
\par
\caption{ 
Simulations (solid lines) versus experiments (symbols) in constant-p triaxial compression
tests on isotropically consolidated samples of Nevada sand \cite{Arulmoli92}.
}%
\label{Fig-Nevada-CIDC}%
\end{figure}
%======================
\begin{figure}[pb]
\begin{center}
{\normalsize 
\mbox{
\subfigure{Experiment} \quad 
\hspace{4cm}
\subfigure{Simulation}
} \\ 
\setcounter{subfigure}{0}%
\mbox{
\subfigure[]{\includegraphics[width=0.48\textwidth]{CIUC-Cyclic-Nevada-test-1.eps}} \quad 
\hspace{-1cm}
\subfigure[]{\includegraphics[width=0.48\textwidth]{CIUC-Cyclic-Nevada-model-1.eps}}
} \\ 
\mbox{
\subfigure[]{\includegraphics[width=0.48\textwidth]{CIUC-Cyclic-Nevada-test-2.eps}} \quad 
\hspace{-1cm}
\subfigure[]{\includegraphics[width=0.48\textwidth]{CIUC-Cyclic-Nevada-model-2.eps}}
}\\
}
\end{center}
\par
\caption{
Simulations versus experiments in undrained cyclic triaxial 
tests on Nevada sand with D$_r\approx$40\% \cite{Arulmoli92}.
 }%
\label{Fig-Nevada-CYUC}%
\end{figure}
%======================
\begin{figure}[pb]
\begin{center}
{\normalsize 
\mbox{
\subfigure[]{\includegraphics[width=0.48\textwidth]{CIDC-Sac-1.eps}} \quad 
\hspace{-1cm}
\subfigure[]{\includegraphics[width=0.48\textwidth]{CIDC-Sac-2.eps}}
} \\ 
}
\end{center}
\par
\caption{
Simulations (solid lines) versus experiments (symbols) in drained triaxial compression
tests on isotropically consolidated  loose and dense samples of Sacramento
river sand \cite{Lee67}.
 }%
\label{Fig-Sacramento-CIDC}%
\end{figure}
%======================


\clearpage


%======================
\begin{figure}[pb]
\begin{center}
\includegraphics[width=0.55\textwidth]{Orlando-Mesh.eps}
\end{center}
\par
\caption{ Illustration of the problem in Example 1 in terms of the soil layering, the finite element mesh, and the input base acceleration. }%
\label{Fig-Isolation}%
\end{figure}
%======================
\begin{figure}[pb]
\begin{center}
{\normalsize \mbox{
\subfigure[Uniform] {\includegraphics[width=0.40\textwidth]{Uniform-StressPath.eps}} \quad
\subfigure[Layered] {\includegraphics[width=0.40\textwidth]{Layered-StressPath.eps}}
}
}
\end{center}
\par
\caption{Variation of shear stress $\sigma_{xz}$ vs vertical effective stress $\sigma_{z}$ for Elements E01 -- E19 during and after shaking.
}%
\label{Fig-Isolation-stresspath}%
\end{figure}
%======================
\begin{figure}[pb]
\begin{center}
{\normalsize \mbox{
\subfigure[Uniform] {\includegraphics[width=0.40\textwidth]{Uniform-StressStrain.eps}} \quad
\subfigure[Layered] {\includegraphics[width=0.40\textwidth]{Layered-StressStrain.eps}}
}
}
\end{center}
\par
\caption{Variation of shear stress $\sigma_{xz}$ vs shear strain $\gamma$ for Elements E01 -- E19 during and after shaking.
}%
\label{Fig-Isolation-stressstrain}%
\end{figure}
%======================
\begin{figure}[pb]
\begin{center}
{\normalsize \mbox{
\subfigure[Uniform] {\includegraphics[width=0.40\textwidth]{Uniform-VoidRatio.eps}} \quad
\subfigure[Layered] {\includegraphics[width=0.40\textwidth]{Layered-VoidRatio.eps}}
}
}
\end{center}
\par
\caption{Variation with time of void ratio $e$ for Elements E01 -- E19 during and after shaking.}%
\label{Fig-Isolation-voidratio}%
\end{figure}
%======================
\begin{figure}[pb]
\begin{center}
{\normalsize \mbox{
\subfigure[Uniform] {\includegraphics[width=0.40\textwidth]{Uniform-au_x.eps}} \quad
\subfigure[Layered] {\includegraphics[width=0.40\textwidth]{Layered-au_x.eps}}
}
}
\end{center}
\par
\caption{Time history of horizontal component of acceleration in solid part of the mixture $a_{u,x}$ for nodes at different elevations during shaking.}%
\label{Fig-Isolation-accel}%
\end{figure}
%======================
\begin{figure}[pb]
\begin{center}
{\normalsize \mbox{
\subfigure[Uniform] {\includegraphics[width=0.48\textwidth]{Uniform-EPP_Contour.eps}} \quad
\subfigure[Layered] {\includegraphics[width=0.48\textwidth]{Layered-EPP_Contour.eps}}
}
}
\end{center}
\par
\caption{Contours of Excess Pore Pressure (kPa) in the soil column at different times.}%
\label{Fig-Isolation-EPPcontour}%
\end{figure}
%======================
\begin{figure}[pb]
\begin{center}
{\normalsize \mbox{
\subfigure[Uniform] {\includegraphics[width=0.48\textwidth]{Uniform-Ru_Contour.eps}} \quad
\subfigure[Layered] {\includegraphics[width=0.48\textwidth]{Layered-Ru_Contour.eps}}
}
}
\end{center}
\par
\caption{Contours of Excess Pore Pressure Ratio ($R_{u}$) in the soil column at different times.}%
\label{Fig-Isolation-Rucontour}%
\end{figure}
%======================


\clearpage


%======================
\begin{figure}[pb]
\begin{center}
\includegraphics[width=0.55\textwidth]{Is-Tokyo-Mesh.eps}
\end{center}
\par
\caption{ Illustration of the problem in Example 2 in terms of the soil layering, the finite element mesh, and the input base acceleration. }%
\label{Fig-mesh-accel}%
\end{figure}
%======================
\begin{figure}[pb]
\begin{center}
{\normalsize \mbox{
\subfigure[]{\includegraphics[ width=0.33\textwidth]{Run-B-01-I-a1-s3-t1-StressPath}} \quad
\hspace{-.5cm}
\subfigure[]{\includegraphics[ width=0.33\textwidth]{Run-B-01-I-a1-s3-t1-StressStrain}}\quad
\hspace{-.5cm}
\subfigure[]{\includegraphics[ width=0.33\textwidth]{Run-B-01-I-a1-s3-t1-Ru}}
} \mbox{
\subfigure[]{\includegraphics[ width=0.33\textwidth]{Run-B-01-I-a1-s3-t1-Ru_Contour}} \quad 
\hspace{-.5cm}
\subfigure[]{\includegraphics[ width=0.33\textwidth]{Run-B-01-I-a1-s3-t1-gamma_Contour}} \quad
\hspace{-.5cm}
\subfigure[]{\includegraphics[ width=0.33\textwidth]{Run-B-01-I-a1-s3-t1-displ_u_x_Contour}} 
}}
\end{center}
\par
\caption{Results of response of the uniform soil column with $a_{max}=0.2g$ (case I) in selected elements along soil profile (E03, E07, E11, E13, E19): (a) shear stress vs. vertical effective stress, (b) shear stress vs. shear strain, (c) excess pore pressure ratio histories, (d) excess pore pressure ratio vs. time and depth, (e) shear strain vs. time and depth, (f) lateral displacement vs. time and depth.}%
\label{Fig-case-I}%
\end{figure}
%======================
\begin{figure}[pb]
\begin{center}
{\normalsize \mbox{
\subfigure[]{\includegraphics[ width=0.33\textwidth]{Run-B-03-I-a1-s3-t1-StressPath}} \quad
\hspace{-.5cm}
\subfigure[]{\includegraphics[ width=0.33\textwidth]{Run-B-03-I-a1-s3-t1-StressStrain}}\quad
\hspace{-.5cm}
\subfigure[]{\includegraphics[ width=0.33\textwidth]{Run-B-03-I-a1-s3-t1-Ru}}
} \mbox{
\subfigure[]{\includegraphics[ width=0.33\textwidth]{Run-B-03-I-a1-s3-t1-Ru_Contour}} \quad 
\hspace{-.5cm}
\subfigure[]{\includegraphics[ width=0.33\textwidth]{Run-B-03-I-a1-s3-t1-gamma_Contour}} \quad
\hspace{-.5cm}
\subfigure[]{\includegraphics[ width=0.33\textwidth]{Run-B-03-I-a1-s3-t1-displ_u_x_Contour}} 
}}
\end{center}
\par
\caption{Results of response of the layered soil column with $a_{max}=0.2g$ (case II) in selected elements along soil profile (E03, E07, E11, E13, E19): (a) shear stress vs. vertical effective stress, (b) shear stress vs. shear strain, (c) excess pore pressure ratio histories, (d) excess pore pressure ratio vs. time and depth, (e) shear strain vs. time and depth, (f) lateral displacement vs. time and depth.}%
\label{Fig-case-II}%
\end{figure}
%======================
\begin{figure}[pb]
\begin{center}
{\normalsize \mbox{
\subfigure[]{\includegraphics[ width=0.33\textwidth]{Run-B-03-I-a2-s3-t1-StressPath}} \quad
\hspace{-.5cm}
\subfigure[]{\includegraphics[ width=0.33\textwidth]{Run-B-03-I-a2-s3-t1-StressStrain}}\quad
\hspace{-.5cm}
\subfigure[]{\includegraphics[ width=0.33\textwidth]{Run-B-03-I-a2-s3-t1-Ru}}
} \mbox{
\subfigure[]{\includegraphics[ width=0.33\textwidth]{Run-B-03-I-a2-s3-t1-Ru_Contour}} \quad 
\hspace{-.5cm}
\subfigure[]{\includegraphics[ width=0.33\textwidth]{Run-B-03-I-a2-s3-t1-gamma_Contour}} \quad
\hspace{-.5cm}
\subfigure[]{\includegraphics[ width=0.33\textwidth]{Run-B-03-I-a2-s3-t1-displ_u_x_Contour}} 
}}
\end{center}
\par
\caption{Results of response of the layered soil column with $a_{max}=0.4g$ (case III) in selected elements along soil profile (E03, E07, E11, E13, E19): (a) shear stress vs. vertical effective stress, (b) shear stress vs. shear strain, (c) excess pore pressure ratio histories, (d) excess pore pressure ratio vs. time and depth, (e) shear strain vs. time and depth, (f) lateral displacement vs. time and depth.}%
\label{Fig-case-III}%
\end{figure}
%======================
\begin{figure}[pb]
\begin{center}
{\normalsize \mbox{
\subfigure[] {\includegraphics[width=0.48\textwidth]{u_x_compare}} \quad
\subfigure[] {\includegraphics[width=0.48\textwidth]{u_x_permanent_compare}}
}
}
\end{center}
\par
\caption{(a) Time histories of surface lateral displacement, and (b) variations with depth of the permanent lateral displacement for cases I, II and III.}%
\label{Fig-ux-compare}%
\end{figure}
%======================


%===============================================
\clearpage

%% The Appendices part is started with the command \appendix;
%% appendix sections are then done as normal sections
%% \appendix

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

\appendix


\section{Constitutive Platform}
\label{Appendix_ConstitutiveModel}


The main components of the SANISAND model are summarized here. Both stress
and strain quantities are assumed positive in compression (as is common in
mechanics of materials), and the effect of this sign convention has been
considered on the model equations. All stress components in this paper
should be considered as effective stress. Finally, in terms of notation, tensor quantities are denoted by bold-faced symbols and operations explained accordingly.  

\textbf{Elasticity:} Isotropic hypoelasticity is defined by%
\begin{equation}
\mathbf{\dot{e}}^{e}=\frac{\mathbf{\dot{s}}}{2G};\ \ \dot{\varepsilon}%
_{v}^{e}=-\frac{\dot{p}}{K}
\label{Eq-Elasticity}%
\end{equation}%
with $\mathbf{s}$\textbf{=}stress deviator; $p$=mean effective pressure; $%
\mathbf{e}$=strain deviator; $\varepsilon _{v}$=volumetric strain;
superscript $e$ denoting elastic; superposed dot denoting the rate; and $G$, 
$K$ the hypoelastic shear and bulk moduli, respectively, given by 
\begin{equation}
G=G_{0}p_{at}\dfrac{(2.97-e)^{2}}{(1+e)}\left( \dfrac{p}{p_{at}}\right)
^{1/2};\ \ K=\dfrac{2(1+\nu )}{3(1-2\nu )}G
\label{Eq-GK}%
\end{equation}%
where $G_{0}$ is a dimensionless material constant, $\nu $ is a constant
Poisson's ratio, $e$ is the void ratio, and $p_{at}$ is the atmospheric
pressure used for normalization.


\textbf{Yield Surface:} The yield surface is defined by 
\begin{equation}
f=[(\mathbf{s}-p\boldsymbol{\alpha }):(\mathbf{s}-p\boldsymbol{\alpha }%
)]^{1/2}-\sqrt{2/3}mp=0  \label{YS}
\end{equation}%
in terms of a deviatoric back-stress ratio $\boldsymbol{\alpha }$, and can
be visualized as a \emph{cone}
The gradient of $f=0$ in stress space is obtained based on Eq. (\ref{YS}) by:%
\begin{equation}
\frac{\partial f}{\partial \boldsymbol{\sigma }}=\mathbf{n}-\frac{1}{2}(%
\mathbf{n}:\mathbf{r})\mathbf{I};\ \ \mathbf{n=}\frac{\mathbf{r-}\boldsymbol{%
\alpha }}{\sqrt{2/3}m}  \label{dfOdsigma}
\end{equation}%
in terms of deviator unit tensor $\mathbf{n}$ ($\mathbf{n:n}=1$) shown along
the radius $\mathbf{r}-\boldsymbol{\alpha }$ of the yield surface circular
trace, where $\mathbf{I}$=the identity
second-order tensor. The $\mathbf{n}$ of Eq. (\ref{dfOdsigma}) is used to
define an effective Lode angle $\theta $ as 
\begin{equation}
\cos {3\theta }=-\sqrt{6}\text{tr}\mathbf{n}^{3}
\end{equation}


The critical stress ratio for a given Lode angle $\theta$ is denoted by $M$ and has been obtained by interpolation between values of $M_{c}$ at $\theta =0$ (compression) and $M_{e}$ at $\theta =\pi /3$ (extension) according to
\begin{equation}
M=M_{c}g(\theta ,c);\ \ g(\theta ,c)=\frac{2c}{(1+c)-(1-c)\cos {3\theta }};\
\ c=\frac{M_{e}}{M_{c}}  \label{DM04Argris}
\end{equation}


\textbf{Bounding, Dilatancy and Critical Surfaces:} Three concentric and
homologous surfaces, the bounding, dilatancy, and critical are considered for the model in the 
$\pi $-plane.
Their tensor points $ \boldsymbol{\alpha }_{\theta }^{b}$, $\boldsymbol{\alpha }_{\theta }^{d}$,
and $\boldsymbol{\alpha }_{\theta }^{c}$ along the direction $\mathbf{n}$
emanating from the origin at a Lode angle $\theta $, are defined
analytically by
\begin{subequations}
\begin{align}
\boldsymbol{\alpha }_{\theta }^{b}& =\sqrt{2/3}[M\exp {(-n^{b}\psi )}-m]%
\mathbf{n}=\sqrt{2/3}\alpha _{\theta }^{b}\mathbf{n} \\
\boldsymbol{\alpha }_{\theta }^{d}& =\sqrt{2/3}[M\exp {(n^{d}\psi )}-m]%
\mathbf{n}=\sqrt{2/3}\alpha _{\theta }^{d}\mathbf{n} \\
\boldsymbol{\alpha }_{\theta }^{c}& =\sqrt{2/3}[M-m]\mathbf{n}=\sqrt{2/3}%
\alpha _{\theta }^{c}\mathbf{n}  \label{alpha_c_theta}
\end{align}
\end{subequations}%
where $\psi =e-e_{c}$=state parameter; $e_{c}=e_{0}-\lambda \left(
p_{c}/p_{at}\right) ^{\xi }$ as the analytical expression of the critical
state line with constants $e_{0}$, $\lambda $ and $\xi $ \cite{Li98}; and $n^{b}$
and $n^{d}$ are material constants.


\textbf{Plastic Flow:} In general the model employs a non-associative flow
rule. The plastic strain rate direction $\mathbf{R}$ will have deviatoric
and volumetric part according to$\mathbf{R}=\mathbf{R}^{\mathbf{\prime }%
}-(1/3)D\mathbf{I}$, where $D$ is the dilatancy coefficient to be specified in the
sequel. The deviatoric part $\mathbf{R}$ is defined as the normal to the
critical surface at point $\boldsymbol{\alpha }_{\theta }^{c}$, Eq. (\ref%
{alpha_c_theta}). The plastic strain is given by 
\begin{subequations}
\begin{align}
\boldsymbol{\dot{\varepsilon}}^{p}& =\left\langle L\right\rangle \mathbf{R}
\label{DM04PF1} \\
\mathbf{R}& =B\mathbf{n}+C(\mathbf{n}^{2}-\frac{1}{3}\mathbf{I})-\frac{1}{3}%
D\mathbf{I} \\
B& =1+\frac{3}{2}\frac{1-c}{c}g\cos {3\theta };\ \ C=3\sqrt{\frac{3}{2}}%
\frac{1-c}{c}g
\end{align}%
where $L$ is the loading index; and $\left\langle {}\right\rangle $ the
Macauley brackets ($\left\langle L\right\rangle =L$ if $L>0$ and $%
\left\langle L\right\rangle =0$ if $L\leq 0$). The dilatancy coefficient $D$
in Equation~(\ref{DM04PF1}) is defined by 
\end{subequations}
\begin{subequations}
\begin{align}
D& =A_{d}(\boldsymbol{\alpha }_{\theta }^{d}-\boldsymbol{\alpha }):\mathbf{n}
\label{DM04PF4} \\
A_{d}& =A_{0}(1+\left\langle \mathbf{z}:\mathbf{n}\right\rangle )
\label{DM04PF5}
\end{align}%
where $A_{0}$ is a material constant. The so-called fabric-dilatancy
internal variable $\mathbf{z}$ has been introduced in order to account for
the effect of fabric changes during the dilatant phase of deformation on the
subsequent contractant response upon load increment reversals. The evolution
law for the fabric-dilatancy tensor $\mathbf{z}$ will be presented later.


\textbf{Evolution Laws:} This model has two tensorial internal variables,
namely, the back-stress ratio tensor $\boldsymbol{\alpha }$ and the
fabric-dilatancy tensor $\mathbf{z}$. The evolution law for the back-stress
ratio $\boldsymbol{\alpha }$ is function of the distance between bounding
and current back-stress ratio in the form of 
\end{subequations}
\begin{equation}
\boldsymbol{\dot{\alpha}}=\left\langle L\right\rangle \frac{2}{3}h(%
\boldsymbol{\alpha }_{\theta }^{b}-\boldsymbol{\alpha })  \label{DM04Evo1}
\end{equation}


While constant values of $h$ yield reasonable simulations, better results
are obtained if $h$ is made function of the current state variables as
\begin{subequations}
\begin{align}
h& =\frac{b_{0}}{(\boldsymbol{\alpha }-\tilde{\boldsymbol{\alpha }}):\mathbf{%
n}}  \label{DM04Evo1a} \\
b_{0}& =G_{0}h_{0}(1-c_{h}e)\left( p/p_{at}\right) ^{-1/2}
\end{align}%
where $\tilde{\boldsymbol{\alpha }}$ is the initial value of $\boldsymbol{\alpha }$ at
initiation of a new loading process and is updated to the new value when the
denominator of Eq.\ (\ref{DM04Evo1a}) becomes negative. The $h_{0}$ and $%
c_{h}$ are material constants.


Finally the evolution law for the fabric-dilatancy $\mathbf{z}$ is
introduced as 
\end{subequations}
\begin{equation}
\mathbf{\dot{z}}=-c_{z}\left\langle -LD\right\rangle (z_{max}\mathbf{n}+%
\mathbf{z})  \label{DM04Evo2}
\end{equation}%
with $c_{z}$ and $z_{max}$ as the material constants that control the
maximum value of $\mathbf{z}$ and its pace of evolution. Eq. (\ref{DM04PF5})
and (\ref{DM04Evo2}) are introduced to enhance the subsequent contraction in
unloading after the dilation in loading in order to generate the well-known butterfly
shape in the undrained cyclic stress path.


%===============================================


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


The overall equilibrium or momentum balance equation for the soil-fluid
\emph{mixture} can be written as 
\begin{equation}
\sigma _{ij,j}-\rho \ddot{u}_{i}-\rho _{f}\ddot{w}_{i}+\rho b_{i}=0
\label{upU_GEequilibrium}
\end{equation}%
where $\ddot{u}_{i}$ is the acceleration of the solid part, $b_{i}$ is the
body force per unit mass, $\ddot{w}_{i}$ is the fluid acceleration relative
to the solid part. For fully saturated porous media (no air trapped inside),
density is equal to $\rho =n\rho _{f}+(1-n)\rho _{s}$, where $n$ is the
porosity, $\rho _{s}$ and $\rho _{f}$ are the soil particle and water
density respectively.


For the pore fluid, the equation of momentum balance can be written as 
\begin{equation}
-p_{,i}-{R_{i}}-\rho _{f}\ddot{w}_{i}-\frac{\rho _{f}\ddot{w}_{i}}{n}+\rho
_{f}b_{i}=0  \label{GEequilibriumFluid}
\end{equation}%
where $R$ is the viscous drag force. According to the Darcy's seepage law,
the viscous drag forces $R$ between soil matrix and pore fluid (water) can
be written as ${R_{i}}=k_{ij}^{-1}\dot{w}_{j}$, where $k_{ij}$ is the tensor
of anisotropic Darcy permeability coefficients. For simple case of isotropic
permeability, scalar value of permeability $k$ is used. The permeability ${k}
$ used here with dimension of [$L^{3}TM^{-1}$] is different from the
permeability used in the usual soil mechanics ($K$) which has the same
dimension of velocity, i.e. [$LT^{-1}$]. Their values are related by $%
k=K/g\rho _{f}$, where $g$ is the gravitational acceleration and the
permeability $K$ is measured in an experiment.


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


In the above governing equations, convective and terms of lower order are
omitted \cite{Zienkiewicz99}. A change of variable is performed by
introducing an alternative variable $U_{i}$, defined as $%
U_{i}=u_{i}+U_{i}^{R}=u_{i}+{w_{i}}/{n}$, that represents absolute
displacement of the pore fluid. The basic set of unknowns is then comprised
of the soil skeleton displacements $u_{i}$, the water pore pressure $p$, and
the water displacements $U_{i}$. These unknown vector ($u_i$ and $U_i$) and
scalar ($p$) fields can be approximated using shape functions and the
corresponding (unknown) nodal values (solid displacements $\bar{u}_{Ki}$,
pore fluid pressures $\bar{p}_K$, and pore fluid displacements $\bar{U}_{Ki}$%
), following standard Finite Element Discretization \cite%
{Zienkiewicz91a,Zienkiewicz91b} 
\begin{equation}
u_i = N_K^u \bar{u}_{Ki}, \;\;\; p = N_K^p \bar{p}_K, \;\;\; U_i = N_K^U 
\bar{U}_{Ki}  \label{39}
\end{equation}
where $N_K^u$, $N_K^p$ and $N_K^U$ are (same) shape functions for solid
displacement, pore pressure and fluid displacement respectively. Each node
of the ($u-p-U$) element has thus seven degrees of freedoms in 3 dimensions
(three for solid displacements, one for pore fluid pressures and three or
pore fluid displacements).
%
After   some   tensor  algebra  and  manipulations, the final $u-p-U$ form can be written as
%
\begin{gather}
\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}^{^{\prime }}\mathrm{d}\Omega  \\ 
0 \\ 
0%
\end{array}%
\right] 
=
\left[ 
\begin{array}{c}
\bar{f}_{Ki}^{u} \\ 
0 \\ 
\bar{f}_{Ki}^{U}%
\end{array}%
\right]  \label{68}
\end{gather}%
The left hand side components of the above matrix equation are given as: 
\begin{eqnarray}
(M_{s})_{KijL}=\int_{\Omega }N_{K}^{u}(1-n)\rho _{s}\delta _{ij}N_{L}^{u}%
\mathrm{d}\Omega \;\;\;\; &\mbox{;}&\;\;\;\;(M_{f})_{KijL}=\int_{\Omega
}N_{K}^{U}n\rho _{f}\delta _{ij}N_{L}^{U}\mathrm{d}\Omega   \nonumber \\
(C_{1})_{KijL}=\int_{\Omega }N_{K}^{u}n^{2}k_{ij}^{-1}N_{L}^{u}\mathrm{d}%
\Omega \;\;\;\; &\mbox{;}&\;\;\;\;(C_{2})_{KijL}=\int_{\Omega
}N_{K}^{u}n^{2}k_{ij}^{-1}N_{L}^{U}\mathrm{d}\Omega   \nonumber \\
(C_{3})_{KijL}=\int_{\Omega }N_{K}^{U}n^{2}k_{ij}^{-1}N_{L}^{U}\mathrm{d}%
\Omega \;\;\;\; &\mbox{;}&\;\;\;\;(G_{1})_{KiM}=\int_{\Omega
}N_{K,i}^{u}(\alpha -n)N_{M}^{p}\mathrm{d}\Omega   \nonumber \\
(G_{2})_{KiM}=\int_{\Omega }nN_{K,i}^{U}N_{M}^{p}\mathrm{d}\Omega \;\;\;\; &%
\mbox{;}&\;\;\;\;P_{NM}=\int_{\Omega }N_{N}^{p}\frac{1}{Q}N_{M}^{p}\mathrm{d}%
\Omega   \label{69-1}
\end{eqnarray}%
while the right hand side components are given as: 
\begin{eqnarray}
(\bar{f}_{s})_{Ki} &=&\int_{\Gamma _{t}}N_{K}^{u}\sigma _{ij}^{^{\prime
}}n_{j}\mathrm{d}\Gamma -\int_{\Gamma _{p}}N_{K}^{u}(\alpha -n)pn_{i}\mathrm{%
d}\Gamma +\int_{\Omega }N_{K}^{u}(1-n)\rho _{s}b_{i}\mathrm{d}\Omega  
\nonumber \\
(\bar{f}_{f})_{Ki} &=&-\int_{\Gamma _{p}}N_{K}^{U}npn_{i}\mathrm{d}\Gamma
+\int_{\Omega }N_{K}^{U}n\rho _{f}b_{i}\mathrm{d}\Omega   \label{69-2}
\end{eqnarray}


%===============================================
\clearpage

%% References
%%
%% Following citation commands can be used in the body text:
%% Usage of \cite is as follows:
%%   \cite{key}         ==>>  [#]
%%   \cite[chap. 2]{key} ==>> [#, chap. 2]
%% 

%% References with bibTeX database:

\section*{References}
\bibliographystyle{elsarticle-num}
%\bibliography{/Users/Mahdi/Library/Bib/Mahdi}
\bibliography{Mahdi}

%% Authors are advised to submit their bibtex database files. They are
%% requested to list a bibtex style file in the manuscript if they do
%% not want to use elsarticle-num.bst.

%% References without bibTeX database:

% \begin{thebibliography}{00}

%% \bibitem must have the following form:
%%   \bibitem{key}...
%%

% \bibitem{}

% \end{thebibliography}

\end{document}

%%
%% End of file `elsarticle-template-num.tex'.
